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Abstract 

Optical  interferometry  is  a  technique  which  can  be  used  to  obtain  high-resolution  imagery 
of  a  distant  target  by  interfering  light  collected  by  multiple  telescopes.  Image  restoration 
from  interferometric  measurements  poses  a  unique  set  of  challenges.  The  first  challenge  is 
that  the  measurement  set  provides  only  a  sparse-sampling  of  the  object's  Fourier  Transform 
and  hence  image  formation  from  these  measurements  is  an  inherently  ill-posed  inverse 
problem.  Secondly,  atmospheric  turbulence  causes  severe  distortion  of  the  phase  of  the 
measured  Fourier  samples.  We  develop  array  design  conditions  for  unique  recovery  of 
the  true  Fourier  phase  in  the  presence  of  this  distortion,  as  well  as  a  comprehensive 
algorithmic  framework  based  on  the  notion  of  redundant-spaced-calibration  (RSC),  which 
together  achieve  reliable  image  reconstruction  in  spite  of  these  challenges.  Within  this 
framework,  we  see  that  the  classical  interferometric  observables  known  as  the  bispectrum 
and  closure  phase  can  limit  sensitivity,  and  that  generalized  notions  of  these  observables 
can  improve  both  theoretical  and  empirical  performance  in  Fourier  phase  estimation.  Our 
framework  leverages  techniques  from  lattice  theory  to  resolve  integer  phase  ambiguities 
in  the  interferometric  phase  measurements,  and  from  graph  theory,  to  select  a  reliable  set 
of  generalized  observables.  As  part  of  the  performance  assessment  of  our  algorithm,  we 
first  show  that  both  the  theoretical  and  simulated  Fourier-phase  estimation  accuracy  of  the 
algorithm  approach  an  atmosphere-oracle  Cramer-Rao  Lower  Bound  at  flux  levels  as  low  as 
60  photons  per  interferometric  fringe  in  the  shot-noise-limited  setting.  Leveraging  techniques 
from  the  field  of  sparse  recovery,  we  then  demonstrate  reliable  image  reconstruction  from  the 
recovered  Fourier  estimates.  Our  results  show  that  reconstructed  image  quality  is  retained 


even  in  simulated  stressing  scenarios  consisting  of  per-exposure  flux  levels  on  the  order  of 
10  photons  per  interferometric  fringe  with  less  than  10  minutes  of  observation  time.  The  end 
result  is  a  comprehensive  strategy  to  achieve  reliable  image  reconstruction  of  dim  objects 
with  optical  interferometry. 
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Introduction 


The  use  of  optical  interferometry  as  a  multi-aperture  imaging  approach  is  attracting  in¬ 
creasing  interest  in  the  astronomical  and  remote-sensing  communities.  The  appeal  of  this 
technique  is  primarily  due  to  the  high  resolution  it  affords  relative  to  single-aperture  imag¬ 
ing.  Namely,  the  diffraction-limited  angular  resolution  of  a  single  aperture  is  proportional  to 
jj,  where  A  is  the  wavelength  of  the  light,  and  D  is  the  diameter  of  the  aperture.  On  the  other 
hand,  the  achievable  angular  resolution  of  an  interferometric  array  of  apertures  is  instead 
proportional  to  g^,  where  Bmax  is  the  maximum  spatial  separation  of  any  two  apertures  in 
the  array.  Therefore  with  interferometry  one  can  achieve  the  same  high  resolution  offered 
by  an  extremely  large  (and  often  prohibitively-costly)  telescope  by  interfering  light  from 
several  telescopes  of  practical  size  distributed  over  a  large  area. 

Image  reconstruction  in  interferometry  is  akin  to  the  more  general  problem  of  imaging 
with  sparse  samples  in  the  Fourier  domain,  which  is  encountered  in  many  other  fields  from 
medical  imaging  to  radio  astronomy.  Namely,  the  fundamental  problem  is  to  regularize 
the  ill-posed  reconstruction  of  NreS£i  resolution  elements  in  the  astronomical  scene  from 
m  <C  NreS£i  Fourier  samples  derived  from  the  array's  interference  patterns.  Traditionally 
approaches  based  on  the  so-called  CLEAN  algorithm  due  to  Hogbom  (1974),  which  implicitly 
fit  the  measurements  to  an  image  model  consisting  of  a  sparse  collection  of  point-like  sources, 
have  been  used  to  solve  this  problem.  Given  the  success  of  this  relatively-simple  algorithm, 
it  is  hardly  surprising  that  recently-developed  techniques  leveraging  sophisticated  sparse 
image  models  (Wiaux  et  al.,  2009)  (Kurien  et  al.,  2014)  have  shown  promise.  In  fact  such 
algorithms  belong  to  a  burgeoning  family  of  cross-disciplinary  techniques,  which  are  based 


1 


on  seminal  work  within  the  last  decade  (Donoho,  2006)  (Candes  et  al.,  2006),  in  the  allied 
fields  of  sparse  recovery  and  compressed  sensing.  These  algorithms  generally  take  advantage 
of  the  fact  that  while  undersampled  Fourier  measurement  sets  do  not  uniquely  define  an 
arbitrary  underlying  image,  they  can  uniquely  specify  a  special  set  of  images  of  practical 
interest:  those  that  are  sparse  in  some  subspace. 

If  Fourier  samples  were  directly  available,  we  could  directly  apply  standard  techniques 
for  image  recovery  from  a  sparse  measurement  set  in  the  frequency  domain.  Flowever, 
estimation  of  the  Fourier  samples  is  itself  a  complicated  inference  problem.  Interferometric 
systems  interfere  signals  from  multiple  apertures  to  produce  superpositions  known  as 
fringes  which  encode  samples  of  the  Fourier  Transform  of  the  object  under  observation.  In 
our  introductory  Chapter  1,  we  present  a  mathematical  model  for  fringe  observation  as 
well  as  sensitivity  limits  for  estimation  of  the  fringe  parameters.  Estimates  of  the  scene's 
brightness-normalized  Fourier  samples,  which  are  known  as  complex  visibilities,  can  then  be 
derived  from  these  parameters.  A  fundamental  challenge  in  interferometry  is  the  distortion 
of  the  phase  of  the  measured  complex  visibilities  due  to  natural  variation  in  the  effective 
path  lengths  to  the  target  observed  by  each  aperture.  In  practice,  the  most  stressing  source 
of  this  variation  is  turbulence  in  the  Earth's  atmosphere.  This  turbulence  alters  the  mean 
phase  profile  at  each  aperture  in  the  array  by  a  non-uniform  and  time-varying  amount  (i.e. 
the  so-called  optical  piston),  which  in  turn  causes  rapid  shifting  of  the  fringes  on  the  focal 
plane.  If  uncorrected,  the  resulting  phase  noise  in  the  Fourier  estimates  then  causes  severe 
distortion  in  the  reconstructed  image. 

Techniques  for  elimination  of  this  phase  noise  fall  into  two  categories  according  to 
the  photon  flux  level.  At  high  photon-flux  levels,  the  fringes  are  bright  and  the  complex 
visibility  estimates  have  high  Signal-to-Noise  ratio  (SNR).  In  this  case,  it  is  possible  to 
directly  decouple  the  contributions  of  the  true  Fourier  phase  from  the  confounding  phase 
contribution  arising  from  the  unknown,  atmosphere-induced  piston  variation.  While  the 
problem  is  ill-posed  in  general,  this  decoupling  can  be  performed  by  alternating  estimation 
of  the  Fourier  and  image  pixels  while  enforcing  a-priori  constraints  in  the  image  domain. 
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such  as  source-sparsity  (Pearson  and  Readhead,  1984).  At  low  flux  levels,  on  the  other  hand, 
the  individual  fringe  exposures  (or  frames )  are  typically  too  weak  for  such  techniques  to 
be  reliable.  Moreover,  the  presence  of  random  phase  variation  across  frames  due  to  the 
atmosphere  means  that  these  fringe  measurements  cannot  be  directly  integrated.  Instead 
atmosphere-invariant  derivatives  1  of  the  fringe  measurements  such  as  the  bispectrum  and 
the  power  spectrum  are  integrated  and  the  complex  visibilities  must  then  be  inferred  from 
these  integrated  observables. 

This  thesis  in  large  part  reflects  the  confluence  of  two  emerging  trends  in  interferometry 
and  in  inverse  imaging  problems  taken  collectively:  in  the  former,  there  is  a  need  for 
approaches  for  imaging  complex,  extended  objects  with  a  sparse  measurement  set  and 
limited  a-priori  knowledge,  and  in  the  latter,  sparse  recovery  and  compressed  sensing  (CS) 
techniques  continue  to  show  great  promise  in  imaging  a  vast  range  of  natural  images  from 
similarly-sparse  measurement  sets.  While  standard  CS  techniques  are  tailored  to  linear- 
inverse  problems  and  do  not  apply  to  atmosphere-invariant  interferometric  observables,  we 
develop  and  validate  an  algorithmic  interface  in  Chapter  2  permitting  these  techniques  to  be 
used  successfully.  We  show  empirically  that  for  compact  scenes  the  effect  of  the  atmosphere 
can  be  eliminated  for  a  wide  range  of  fluxes.  The  limiting  restriction  of  our  interface  to  the 
imaging  of  compact  scenes  is  a  direct  consequence  of  the  fundamental  ill-posed  nature  of 
Fourier  phase  recovery  from  atmosphere-perturbed  fringe  measurements.  Namely,  since 
the  mapping  from  Fourier  phase  to  the  set  of  measured  phases  is  rendered  non-injective  2 
by  the  atmosphere  in  general,  we  must  introduce  constraints  upon  the  object  to  uniquely 
determine  the  Fourier  phase.  Such  regularization  techniques  and  the  object  constraints  they 
impose  can  be  obviated  by  introducing  redundancy  into  the  inter-aperture  spacings  in  the 
interferometric  array.  This  well-posed  framework  for  Fourier  phase  determination  forms  the 
basis  of  Chapters  3  and  4. 

xBy  derivative,  we  mean  an  observable  derived  from  the  fringe  measurements,  and  not  the  mathematical 
derivative. 

2 By  non-injective,  we  simply  mean  that  a  given  measurement  set  can  be  produced  by  different  sets  of  Fourier 
phases  so  that  the  former  does  not  uniquely  determine  the  latter. 
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The  use  of  redundant  aperture  spacing  to  recover  the  Fourier  phase  in  the  presence  of 
piston  variation  in  a  well-posed  manner  is  known  as  redundant  spacing  calibration  (RSC).  Of 
the  myriad  RSC  techniques  that  have  been  developed  for  both  radio  and  optical  interferom¬ 
etry,  some  operate  on  the  measured  complex  visibilities  (which  we  call  Phasor  approaches), 
and  others  on  the  phase  component  of  these  complex  visibilities  (which  we  call  Phase 
approaches).  Even  with  redundancy,  a  fundamental  ambiguity  exists  in  RSC-based  phase 
recovery  which  is  rooted  in  the  27r-periodicity  of  the  interferometric  phase.  As  a  direct  result 
of  this  ambiguity,  the  solution  for  the  complex  visibilities  given  the  measured  visibilities  is 
non-unique  in  both  Phase  and  Phasor  approaches.  Flowever,  we  will  show  in  Chapter  3  that 
for  certain  patterns,  which  we  denote  wrap-invariant  patterns,  this  fundamental  ambiguity 
can  be  rendered  benign.  Namely,  for  such  patterns,  phase-unwrapping  techniques  such  as 
those  based  on  the  closest-vector-problem  (CVP)  formulation  due  to  Lannes  and  Anterrieu 
(1999)  can  successfully  recover  the  true  Fourier  phase  (modulo  2n).  We  show  that  this 
wrap-invariance  property  is  conferred  upon  arrays  whose  interferometric  graph  satisfies 
a  certain  cycle-free  condition.  This  condition  is,  to  the  best  of  our  knowledge,  the  first 
sufficient  condition  on  an  interferometric  aperture  pattern  for  unique  recovery  of  the  Fourier 
phase.  For  cases  in  which  this  condition  is  not  satisfied,  we  provide  a  simple  algorithm 
for  identifying  those  graph  cycles  which  prevent  its  satisfaction.  For  illustrative  purposes, 
we  apply  this  algorithm  to  diagnose  a  member  of  a  aperture-pattern  family  popular  in  the 
literature  which  is  not  wrap-invariant,  and  modify  it  so  that  it  achieves  wrap-invariance. 

Flaving  established  conditions  for  unique  Fourier  phase  recovery,  we  then  turn  our 
attention  in  Chapter  4  to  the  issue  of  well-posed  and  practical  image  formation.  In  partic¬ 
ular,  we  consider  low-flux  scenarios  in  which  atmosphere-invariant  observables  must  be 
integrated  over  many  frames  for  reliable  phase  estimation.  In  hopes  of  improving  sensitivity, 
we  generalize  the  classical  atmosphere-invariant  phase  observables  (i.e.  the  bispectrum 
and  closure  phase)  to  higher-order  observables,  which  we  denote  as  the  n-spectrum  and 
generalized  closure  phase  respectively.  We  extend  the  uniqueness  results  in  Chapter  3  to  our 
new  observables,  and  develop  a  novel  comprehensive  algorithm  for  image  reconstruction 
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from  these  observables.  Here  we  leverage  the  notion  of  minimum  cycle  basis  from  graph 
theory.  A  standard  sparse-recovery  technique  known  as  total-variation  minimization  is  used  to 
perform  the  final  image  reconstructions.  We  present  numerical  and  visual  evidence  indicat¬ 
ing  that  our  generalized  observables  can  yield  better  estimation  performance  relative  to  their 
classical  counterparts.  Moreover,  we  show  that  our  algorithm's  performance  approaches  the 
Cramer-Rao  Lower  Bound  for  this  estimation  problem  in  an  example  scenario  of  imaging  a 
dim  object. 

In  summary,  great  progress  has  been  made  recently  in  solving  standard,  linear  inverse 
problems  in  imaging  as  part  of  the  development  of  compressed  sensing  and  sparse  recovery 
methodology.  Simultaneously,  many  algorithmic  techniques  have  been  developed  for 
solving  the  specialized  ill-posed  inverse  problem  in  optical  interferometry,  which  is  greatly 
complicated  by  the  atmosphere.  Characterizing  the  penalty  in  achievable  performance 
associated  with  this  increased  complexity  has  remained  an  open  problem.  Hence  in  this 
thesis  we  ask  the  question:  "To  what  extent  can  we  make  interferometric  imaging  immune  to 
the  effects  of  the  atmosphere  in  both  reliability  and  sensitivity?".  We  then  probe  the  limits  of 
this  immunity  using  a  novel  algorithmic  framework  which  guarantees  unique  phase  recovery 
through  the  notion  of  wrap-invariance,  and  generalizes  the  notion  of  atmosphere-invariant 
observable  in  order  to  approach  theoretical  sensitivity  limits  at  low  flux  levels.  In  these  two 
pursuits,  our  framework  leverages  theory  and  algorithms  from  lattice  theory  and  graph 
theory,  respectively.  In  the  end  we  provide  strong  evidence  that  even  in  scenarios  thought 
to  be  stressing  from  an  SNR  perspective,  we  can  achieve  near-immunity  to  the  effects  of  the 
atmosphere. 
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Chapter  1 


Fundamentals  of  Optical 
Interferometry 

1.1  Chapter  Overview 

In  this  chapter,  we  introduce  the  physics-based  principles  of  optical  interferometry,  thereby 
providing  a  foundation  for  the  main  results  of  this  thesis.  As  the  cornerstone  of  these 
principles,  the  Van-Cittert-Zernike  Theorem  establishes  the  interferometer  as  a  Fourier- 
imaging  device,  i.e.  one  whose  output  can  be  easily-conceptualized  when  examined  in  the 
Fourier  domain.  In  particular,  we  show  that  interferometric  arrays  characterize  the  object 
under  observation  via  measurement  of  its  brightness-normalized  Fourier  components,  which 
are  the  so-called  complex  visibilities.  We  then  provide  a  derivation  of  the  sensitivity  limits  for 
measurement  of  these  complex  visibilities  in  the  shot-noise-limited  case.  This  derivation 
is  based  on  the  well-known  Cramer-Rao  Lower  Bound  (CRLB)  in  estimation  theory,  and 
corroborates  previous  results  due  to  Zmuidzinas  (2003).  This  CRLB  will  prove  useful  as  a 
performance  benchmark  for  the  algorithms  developed  in  this  thesis.  We  then  proceed  to 
introduce  the  role  of  atmospheric  turbulence  and  other  sources  of  interferometric  phase 
noise  in  interferometry.  Finally  we  define  a  few  important  parameters  of  an  interferometer: 
resolution  element,  Field-of-Viezv,  and  undersampling  ratio. 
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1.2  Interferometric  Architectures 


The  main  appeal  of  interferometry  over  other  astronomical  imaging  techniques  is  that 
it  affords  the  same  high-resolution  capability  offered  by  an  extremely  large  (and  often 
prohibitively-costly)  telescope  by  interfering  light  from  several  telescopes  of  practical  size. 
Optical  interferometers  enable  the  imaging  of  a  distant  scene  by  providing  a  sampling  of 
the  scene's  the  2D  Fourier  Transform.  Several  excellent  surveys  provide  the  theoretical 
basis  of  interferometry  as  well  as  the  practical  issues  involved  in  building  and  operating  an 
interferometer,  including  those  by  Labeyrie  et  al.  (2006),  Glindemann  (2011),  and  Buscher 
(2015).  Each  pair  of  telescopes  in  an  interferometric  array  measures  a  single  angular  spatial 
frequency  of  2^b  radians,  where  b  is  the  vector  difference  of  the  telescope  positions,  which 
is  known  as  a  baseline.  For  an  array  of  Nap  apertures,  the  data  set  then  consists  of  all  (A^p) 
such  measurements. 

There  are  two  popular  beam  combination  architectures  in  use  in  optical  interferometry: 
the  pairwise  combination  scheme,  and  the  Fizeau  combination  scheme.  The  two  schemes 
are  illustrated  in  Figure  1.1.  Suppose  we  have  a  telescope  array  consisting  of  Nap  apertures, 
each  of  which  collects  n  photons  from  a  distant  source.  In  the  Fizeau  scheme,  light  from 
all  telescopes  is  interfered  on  a  single  focal  plane,  forming  an  interference  pattern  known 
as  a  fringe  for  each  telescope  pair.  Each  fringe  encodes  a  sample  of  the  scene's  2D  Fourier 
Transform.  In  the  pairwise  scheme,  light  from  each  aperture  is  split  Nap  ways  and  combined 
with  that  of  each  other  aperture  to  form  a  single  fringe  pattern  on  separate  focal  planes. 
Hence  each  of  the  (N“f’)  focal  planes  receives  N2n_1  a  photons,  which  is  typically  a  small 
fraction  of  the  total  photons  incident  upon  the  array. 

Ideally  an  interferometer  1  provides  a  perfect  encoding  of  the  scene's  sampled  Fourier 
Transform.  In  practice,  however,  we  never  observe  pristine  interference  fringes  in  either 
architecture.  The  quality  of  the  fringes  is  degraded  by  statistical  fluctuations  in  the  arrival 
times  of  photons  on  the  focal  plane,  which  is  a  phenomenon  known  as  shot  noise.  We  will 

1We  will  use  the  terms  interferometer  and  interferometric  array  interchangeably  in  this  thesis 
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Figure  1.1:  The  two  popular  beam  combination  schemes  in  optical  interferometry 


quantify  the  impact  of  shot  noise  further  in  the  course  of  the  thesis.  For  now,  we  note 
that  the  Fizeau  architecture  uses  all  light  to  form  each  fringe  while  incurring  shot  noise 
from  all  N„pn  photons.  On  the  other  hand,  the  pairwise  architecture  uses  a  fraction  of  the 
light  to  form  each  fringe,  which  is  then  by  corrupted  by  the  shot  noise  due  to  the  photons 
collected  by  two  apertures.  A  fundamental  result  due  to  Zmuidzinas  (2003)  quantified 
this  tradeoff  and  established  the  superior  overall  sensitivity  of  the  Fizeau  (or  all-in-one ) 
scheme  with  respect  to  the  pairwise  scheme.  Because  the  Fizeau  architecture  is  also  typically 
simpler  to  implement,  it  is  often  preferred  in  practice.  The  results  in  this  thesis  are  therefore 
geared  toward  the  Fizeau  architecture,  although  we  will  leverage  the  pairwise  architecture 
in  Chapter  4  as  a  conceptual  springboard  for  subsequent  analysis. 


1.3  The  Van-Cittert-Zernike  Theorem 

To  derive  the  Theorem,  let  us  model  a  source  in  the  sky  as  a  superposition  of  differential 
emitting  elements  of  size  AO.  Consider  then  the  electric  field  at  the  y'-th  aperture  of  an 
interferometric  array  which  is  a  distance  d  from  one  such  differential  patch.  This  can  be 
written  as: 


£;(  fi) 


_  fl(0)AO  Juj(t — L)\ 

R:  C  > 


(1.1) 
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Figure  1.2:  The  Fizeau  Interferometer  Concept 


where  a  (Cl)  is  the  amplitude  of  the  field  at  the  source,  Rj  is  the  distance  from  the  object 
to  the  j-th  aperture,  to  is  the  angular  frequency  of  the  wave,  and  c  is  the  speed  of  light.  For 
the  purposes  of  this  analysis,  we  neglect  the  amplitude  variation  amongst  elements  of  the 
array,  so  that  the  amplitude  at  all  apertures  is  a  constant  a(njAn. 

The  fields  from  all  of  the  point  sources  in  the  interferometer's  field-of-view  are  super¬ 
imposed  on  the  system's  focal  plane.  For  the  Fizeau  interferometer,  this  superposition  is 
typically  performed  by  a  beam-combiner,  which  focuses  the  light  collected  by  each  aperture 
to  a  common  point  with  propagation  direction  specified  by  an  aperture-specific  wave-vector 
ky  (see  Figure  1.2). 

Using  the  definitions  above,  we  can  compute  the  field  contribution  of  each  aperture  at 
an  arbitrary  point  on  the  focal  plane  p.  This  contribution  will  be  a  function  of  the  respective 
wave-vectors  and  vector  distances  from  p.  Namely,  we  have: 

Ej(p)  =  a(QjA°  •  ei(ky*,-(p))e^( t-%)  (1.2) 

where  k;  is  the  wavevector  of  the  outgoing  wave  from  aperture  j  at  the  beam-combiner, 
and  Xj(p)  is  the  vector  position  to  p  relative  the  j-th  aperture  position  at  the  beam-combiner. 
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As  we  show  in  Appendix  A.l,  we  can  rewrite  this  expression  in  terms  of  the  lateral 


aperture  position  ry  at  the  beam  combiner  as: 


„  ,  ,  fl(Q)AQ 

EM  =  -4— 


Krj-p+<Pj)eM  *--£). 


■  eK  ’ 1 


(1.3) 


where  c pj  is  an  aperture-dependent  phase  which  is  independent  of  p. 

Now  consider  the  field  superposition  at  the  pixel  at  vector  coordinate  p  on  the  focal 
plane  (see  Figure  1.2),  which  can  be  written  as: 


E( Q)  =  .  (£jV  g*'( rrP)eMt-^))  (1.4) 

where  we  have  neglected  any  attenuation  incurred  by  the  field  en-route  to  the  focal  plane  for 
simplicity.  To  obtain  the  field  for  the  complete  source,  we  now  integrate  over  all  differential 
patches  (which  we  call  point  sources )  in  the  source  to  obtain: 

Etotai  =|dfi-^(f;  (1.5) 

The  detector  at  pixel  p  measures  the  intensity  Q  of  the  superposition.  Namely, 

Q(p)  =  \Efotal\  =  E-totalE-total  (1-6) 


1  p  R-  p  L^aV  r 

Q(p)  =  -^  dClt-  a{ CL1)(Y,eiu(t~^)ei(r>p+,pi))  /  dCt2  •  a*(D2)(£  e-«'«(f-4)e-ife-P+^)) 

®  ■'  7=1  ■’  k=  1 

The  final  measurement  is  a  photon  count  which  is  proportional  to  the  time-averaged 
intensity.  Hence,  after  re-arranging  and  time-averaging,  we  obtain: 


1  p  p  R-  L^aP  R 

( q(p ))  =  -j2  / /  (fl(ni)fl*(n2))(^eia,(f“-: V(r;-/»+^))(^g-i«(t-^.)g-i(*i-p+fc))dn1dn2 

“  dJ  y=i  fc=l 

(1.7) 
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Figure  1.3:  Path  difference  between  two  apertures 


We  can  simplify  this  expression  greatly  by  noting  that,  with  a  few  exceptions,  astronomi¬ 


cal  sources  are  spatially-incoherent,  which  means  that  the  complex  amplitudes  of  radiated 
waves  from  different  points  are  uncorrelated.  Namely, 


(fl(Qi)a*(Q2))  =  -  02) 


(1.8) 


where  d(x)  denotes  the  Dirac  delta  function. 

The  result  is  that  the  double  integral  above  collapses  to  a  single  integral  in  Q.  We 
can  now  write  the  sum  of  all  pairwise  products  in  Equation  (1.7)  as  the  following  single 


summation  over  all  ( N“p )  aperture  pairs: 


where  the  first  term  represents  the  sum  of  the  products  of  each  aperture's  field  with  its 
conjugate  (i.e.  terms  of  the  form  EjEj). 

To  understand  the  meaning  of  Equation  (1.9),  it  is  useful  to  approximate  A using  the 
geometry  illustrated  in  Figure  1.3.  Let  us  define  the  unit  vector  originating  at  aperture  1 
and  pointing  in  the  direction  of  the  emitting  source  as  n.  For  an  arbitrary  source  location. 
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this  unit  vector  would  of  course  differ  for  each  aperture.  In  most  astronomical  scenarios, 
however,  the  inter-aperture  distance  is  extremely  small  compared  with  the  distance  between 
the  source  and  the  aperture  array  (Buscher,  2015).  In  such  cases  it  is  reasonable  to  assume 
that  the  direction  vector  n  is  approximately  constant  across  the  array.  As  shown  in  Figure 
1.3,  we  have: 


AKj*  =  b jk  •  n  (1.10) 

where  b  ^  is  the  vector  difference  position  of  the  k- th  and  j-th  apertures. 

Moreover,  we  can  write  the  differential  patch  dCl  in  terms  of  components  of  this  vector 
n  as: 

dd  =  ( d2)dnxdny  (1.11) 

Substituting  Equations  (1.10)  and  (1.11)  into  Equation  (1.9)  and  noting  that  7  =  X'  we 
obtain: 


(Q(p))  —  J  I(Cl)Napdnxdnv+ 

Ap  />  27Tb  jk  r  27rb jk 

£  (e^-Ar;*+^)  I(Ct)el^h)dnxdny  +  e-l^^+^  l{0)e~l^^dnxdny)  (1.12) 

(j,k):k>j 

Note  that  the  two  integrals  are  F*(x)  and  F(x),  respectively,  where  F(x)  the  2D 

b 

Fourier  Transform  of  the  source  evaluated  at  spatial  frequency  -j-.  For  simplicity  of  notation, 
let  us  index  all  of  the  aperture  pairs  (j,  k)  with  a  single  index  h.  Defining  (pi,  :=  (pjk,  and 
Fj,  :  =  F(  ^ )  and  substituting,  we  have: 


.  (N2P) 

(Q(p))  =  /  I(Cl)Napdnxdny  +  £  Fhe^  Al< *+^)  +  (1.13) 

J  h= 1 

(Q(p))  =  F0Nap  +  £  2|F;,|cos(p  •  Arh  +  ZFh  +  <j>h)  (1.14) 

h= 1 

From  this  last  equation,  we  see  that  focal  plane  consists  of  a  superposition  of  2D  sinusoids 
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(or  fringes),  each  of  which  encodes  a  distinct  Fourier  component  of  the  source.  This  result  is 
the  Van-Cittert-Zernike  Theorem  which  was  first  derived  by  Pieter  van  Cittert  in  1934  (van 
Cittert,  1934)  and  then  proved  in  simpler  fashion  by  Frits  Zernike  in  1938  (Zernike,  1938). 

For  subsequent  analysis  in  this  thesis,  it  will  be  convenient  to  parameterize  the  fringes 
using  the  following  Definitions: 

Definition  1.3.1.  The  complex  fringe  phasor  z  at  a  given  spatial  frequency  is  given  by:  z  := 

\Fb\e’ZFb 

Note  that  the  magnitude  of  the  fringe  phasor  will  be  proportional  to  the  brightness  of 
the  object.  It  is  also  useful  to  have  a  quantity  describing  the  strength  of  a  given  Fourier 
component  relative  to  the  overall  brightness  of  the  object.  This  relative  measure  is  provided 
by  the  complex  visibility: 

Definition  1.3.2.  The  complex  visibility  v  of  an  object  at  a  given  spatial  frequency  is  given  by 

,1  ,•  |F;,|dZFfc 

the  ratio  v  =  J-W — 

to 

Definition  1.3.3.  The  visibility  7  of  an  object  at  a  given  spatial  frequency  is  the  modulus  of 
the  complex  visibility,  i.e.  7  :=  Ff 

Definition  1.3.4.  The  Fourier  phase  6  of  an  object  at  a  given  spatial  frequency  is  the  argument 
(i.e.  phase)  of  the  complex  visibility,  i.e.  9  :=  ZF]V 

With  Definition  1.3.3,  we  can  write  Equation  (1.15)  as: 

(N“p\ 

{ Q(P ))  =  FoKv  +  27 hF0cos(p  ■  Arh  +  ZFh  +  <ph)  (1.15) 

h= 1 

The  expected  photon  counts  observed  at  the  detectors  are  directly  proportional  to  the 
intensity  present.  Let  the  proportionality  constant  relating  intensity  to  photon  counts  be 
given  by  k,  so  that  a  vectorized  representation  of  these  expected  photon  counts  is  given  by: 

(Nap) 

(y (p))  =  KF0Nap  +  k  27 hF0cos(p  ■  Arh  +  ZFh  +  (ph)  (1.16) 

/!  =  1 
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Moreover  we  assume  that  light  is  conserved  throughout  the  fringe  generation  process, 
and  hence  the  photon  counts  must  sum  to  the  total  number  T;,  of  photons  incident  upon  the 
array,  i.e.  Tp  =  nNapr  i.e. 


T> 


E<y(p)>  =  E 

p  p 


/  (N?)  \ 

KF0Nap  +  K  E  27hFocos(p  ■  Ar,;  +  ZFh  +  <ph) 

\  h=1  J 


nNap  (1.17) 


In  this  thesis,  we  will  consider  the  scenario  in  which  the  periods  of  the  fringe  sinusoids 
are  all  integer  multiples  of  pixels  on  the  focal  plane.  This  is  one  of  the  so-called  DFT 
conditions  commonly  employed  in  interferometry  to  avoid  fringe  estimation  bias  due  to 
spectral  leakage  (Gordon,  J.  A.  and  Buscher,  D.  F.,  2012).  In  this  case  the  sinusoidal  term 
vanishes  in  the  sum  in  Equation  (1.17)  leaving  only  the  constant  (first)  term,  and  if  we  define 
Nfp  as  the  number  of  pixels  on  one  side  of  a  square  focal  plane,  we  have: 


Y^KFoNap  =  NfpKFoNap  =  nNap  (1.18) 

p 

which  in  turn  implies  that  K  =  F  ”  2  ,  and  hence  by  substitution  into  Equation  (1.16) 
yields: 

(Naf>\ 

nN  2 n  y  2 

(y(p)>  =  -T#  +  7T  £  7 hCOs(p  •  Ar;,  +  eh  +  <ph)  (1.19) 

J  7p  h= i 

1.4  The  Cramer-Rao  Bound  for  the  Complex  Visibilities 

There  are  two  principal  sources  of  noise  that  affect  this  measurement.  We  have  already 
alluded  to  the  shot  noise  arising  from  the  fact  that  photon  arrivals  are  not  uniform  in  time 
even  if  incident  intensity  is  constant.  The  resulting  uncertainty  in  photon  counts  for  a 
given  exposure  time  is  well-modeled  as  a  Poisson  distribution  with  mean  AT,  where  A  is 
proportional  to  the  intensity  and  T  is  the  exposure  time.  The  second  source  of  noise  is  read 
noise,  which  is  thermal  noise  in  the  detection  circuitry.  Read  noise  is  well-modeled  by  the 
Gaussian  distribution.  Since  shot  noise  is  proportional  to  intensity  whereas  read  noise  is 
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not,  the  former  dominates  the  latter  in  high  light-level  scenarios.  And  with  recent  advances 
in  electronics,  the  light-level  regime  in  which  read  noise  is  important  continues  to  diminish. 
We  will  therefore  focus  on  detection  limits  for  the  shot-noise-limited  case. 

Recall  Equation  (1.19)  provides  an  expression  for  the  time-averaged  intensity.  To  obtain 
lower  bounds  on  sensitivity,  let  us  assume  that  the  phase  differences  A 0;,  are  known  and 
therefore  calibrated.  We  can  write  the  observation  model  in  matrix  form  as: 


W-Nap  2il 


(yW>=^r^ 


cos(p(l)  •  Aki) 
cos(p( 2)  •  Aki) 


— sin(p(l)  •  Aki) 
—sin(p(2)  ■  Aki) 


cos(p(  1)  •  Ak2) 
cos(p(2)  ■  Ak2) 


— sm(p(l)  •  Ak2) 
—sin(p( 2)  •  Ak2) 


Re[v  i] 
Im[vi\ 
Re[v2 ] 
Im[v2\ 


Suppose  that  the  interferometer's  output  is  a  beam  with  an  extent  of  Nfp  pixels  on 
the  focal  plane.  Letting  A  be  the  Nj  p-by-2(N^p)  matrix  in  the  equation  above  (the  so-called 
visibility -to-pixel  matrix  (V2PM)),  we  can  write  this  matrix  equation  compactly  as: 


y 


nNaf 

~*Tr 


■1 


2n  *  - 

NfP  +  ^fpAv 


(1.20) 


where  denote  the  constant  ones  vector  of  length  Nfp,  v  is  the  vector  containing  the 
quadrature  components  of  the  complex  visibilities.  The  actual  pixel  counts  can  be  modeled 
as  i.i.d.  Poisson  random  variables  with  means  y,  i.e. 


Y  ~  Poisson(y).  (1.21) 

Following  a  similar  derivation  given  in  Harmany  et  al.  (2012),  we  can  now  write: 
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(1.22) 


Nfp  'Yi 


p(Yly)  =  n  exv(-eJy) 


i= 1 


Y/! 


where  e,  is  the  unit  vector  of  the  z-th  canonical  basis.  Substituting  for  y,  we  obtain  the 
negative  log-likelihood  as: 


N. 


'fp 


F (y)  =  lTNfpy  -  E  Yil°g(.eI y)  +  c 


(1.23) 


!=1 


where  1  is  an  all-ones  vector  of  size  Nfp,  and  C  is  a  constant  independent  of  v.  Let 


A  :=  ^f-A.  Then  the  gradient  with  respect  to  v  is  given  by: 


NfP 


N.fp  y. 

VvF(y)  =  AT1  -  £  -y-Are,  (1.24) 

Therefore  the  Hessian  is: 

Nfp  y. 

V|F(y)  =  Ar  E  7^2eieiA  (L25) 

i=i  (c  y) 

To  obtain  the  Fisher  information  matrix  I(v),  we  take  the  negative-expectation  of  this 
quantity: 


Nf; 


I(v)  =  -E[V|F(y)]  =  -Ar  f  ^2^ A 

i= i  (e/y) 


(1.26) 


But  since  E[Y,]  =  ejy,  this  simplifies  to: 


I(v)  =  AtDA 


(1.27) 


where  D  is  a  diagonal  matrix  with  D (!  =  prpy- 

Therefore  the  Cramer  Rao  Lower  Bound  (CRLB)  for  the  variance  of  each  estimated 
Fourier  coefficient  is  given  by: 
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Figure  1.4:  The  relationship  between  aperture  patterns  and  Fourier  sampling 


var(vi)  >  [I(v)  X]u  (1.28) 

This  expression  matches  the  result  given  by  Zmuidzinas  (2003). 

1.5  Key  Parameters  of  an  Interferometer:  Resolution,  Field-of-View, 
and  Undersampling  Ratio 

In  the  Section  1.3  we  established  that  the  interference  fringe  generated  by  a  pair  of  apertures 
separated  by  a  vector  b  encodes  the  amplitude  and  phase  of  the  Fourier  Transform  of 
the  scene  at  spatial  frequency  For  real  images  each  baseline  actually  contributes  two 
spatial  frequency  samples,  as  the  Fourier  Transform  of  real  images  is  conjugate  symmetric. 
As  an  example,  consider  the  aperture  pattern  in  the  left  panel  of  Figure  1.4.  Suppose  we 
are  interfering  light  at  a  wavelength  of  500  nm.  The  baseline  bi2  formed  by  apertures 
1  and  2  generates  a  fringe  encoding  the  complex  visibility  V12  at  a  spatial  frequency  of 
=  (2<?6, 6<?6)  cycles /radian,  and  its  conjugate  v(2  at  a  spatial  frequency  —  (2e6, 6e6) 
cycles/ radian. 

Consider  the  Cartesian  frequency  sampling  grid  with  gridpoint  spacing  A as  shown 
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Figure  1.5:  Representation  of  a  scene  as  a  Field-of-Vieiv  comprised  of  resolution  elements 


in  Figure.  Suppose  we  sample  a  2D-signal  in  the  Fourier  domain  with  sampling  interval 
A s,freq  =  jp  along  each  dimension.  By  the  well-known  Nyquist  Sampling  Theorem,  we  know 
that  only  images  spatially-limited  to  ±R  along  each  dimension  can  be  uniquely  represented 
with  these  samples;  images  beyond  this  extent  will  suffer  from  aliasing.  Therefore  an 
interferometer  whose  sampled  spatial  frequencies  lie  on  a  grid  with  spacing  A =  l-ff- 
has  an  unambiguous  Field-of-View  area  of  }- -  =  .  Conversely,  by  the  space-frequency 

° min  umin 

duality  of  the  Nyquist  Theorem,  we  know  that  we  can  uniquely  represent  an  image  ban- 
dlimited  to  frequency  extent  ±L  by  sampling  in  space  at  an  interval  of  As, space  =  4p-  Given 
the  maximum  spatial  frequency  observed  by  our  interferometer  is  L  =  Amax  =  the 
bandlimited  approximation  of  the  image  is  uniquely  defined  by  samples  spaced  jpfp  apart 
along  each  dimension.  Therefore  the  size  of  the  smallest  resolvable  element  in  the  image, 
which  we  call  the  resolution  element,  is  given  by  .  Figure  1.5  depicts  the  notions  of 
Field-of-View  (FOV)  and  resolution  element  (resEl)  visually. 

With  the  areas  of  the  Field-of-View  and  resolution  element  in  hand,  we  can  now  compute 
the  number  of  resolution  elements  in  the  bandlimited  approximation  of  the  image  as  the 
ratio  of  these  two  quantities,  or: 


NresEls  — 


(1.29) 
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where  r  =  ^M2L. 

umin 


We  now  define  our  undersampling  p  as  the  ratio  of  the  number  of  distinct  Fourier 
samples  to  the  number  of  resolution  elements  NreS£is  in  the  image.  Each  Fourier  sample 
we  measure  is  actually  a  pair  of  samples  since  for  real  images,  the  value  of  the  Fourier 
Transform  at  a  given  spatial  frequency  is  the  conjugate  of  that  at  the  corresponding  negative 
spatial  frequency.  Hence  we  have: 


_  2(Nip )  _  1  (N„p\ 

^  NresEls  2r2  \  2  J 


(1.30) 


1.6  Atmospheric  and  Instrumental  Phase  Noise 

Our  analysis  thus  far  has  assumed  that  the  effective  path  between  the  beam-combiner  and  the 
target  depends  merely  on  the  source-instrument  geometry  depicted  in  Figure  1.2.  In  practice 
there  are  several  deterministic  and  random  sources  of  path  variation  across  the  aperture  array. 
Common  deterministic  sources  of  path  variation  include  that  in  the  optical  paths  traversed 
between  the  instrument's  front-end  and  beam-combiner.  The  predominant  random  sources 
of  delay  are  atmospheric  turbulence  and  instrumental  vibration.  Atmospheric  turbulence, 
which  is  typically  the  more  stressing  of  these  two  sources,  alters  the  optical  path  traversed 
by  the  wave  arriving  at  each  aperture  in  a  nonuniform  and  time-varying  manner. 

Let  us  consider  the  total  path  alteration  due  to  all  of  the  sources  mentioned.  Let  us  then 
define  the  corresponding  phase  shift  resulting  from  this  alteration  at  aperture  j  as  el^Mal .  If 
these  phase  shifts  are  carried  through  the  analysis  in  Section  1.3,  we  obtain: 

(N„p) 

nN  2 n  y  2 

(y (P))  =  -i-ff-  +  ttt  E  7 jkcos(p  ■  Arjk  +  0jk  +  (pjk)  (1.31) 

J  Vp  iV/p  m 

where  c pjk  :  =  (pjrt0tal  ~  <pk, total  is  the  difference  between  the  total  phase  shifts  at  the  two 
apertures  in  the  baseline  associated  with  aperture  pair  (/,/). 

Let  us  now  examine  a  traditional  method  for  eliminating  the  effect  of  this  phase  noise 
which  was  first  suggested  by  Jennison  (1958)  in  the  context  of  inteferometry  at  radio  wave- 
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lengths.  Suppose  the  atmosphere  adds  phases  < pj  and  fa  to  apertures  k  and  j,  respectively. 
The  result  is  that  the  observed  complex  visibility  is  given  by: 

vyfc  =  Vfl'jtorM  (1.32) 

We  can  eliminate  such  nuisance  factors  by  forming  a  triple  product  g  of  the  Fourier 
phasors  along  a  triangle  of  baselines  (e.g.  gi23  :=  viiWsWi)-  As  shown  in  Figure  1.6, 
this  special  triple  product  (known  as  the  bispectrum)  cancels  the  atmospheric  phase  terms, 
leaving  only  the  desired  Fourier  information  (i.e.  the  {0}).  The  phase  of  the  bispectrum 
Zg,  which  is  known  as  the  closure  phase,  can  hence  be  used  to  recover  estimates  for  the 
Fourier  phases.  Consider  an  interferometer  that  measures  all  baselines  among  Nap  apertures. 
Out  of  the  (NjP)  possible  triangles,  only  ( 1 )  of  the  associated  closure  phase  relations 
are  linear ly-independent  (Readhead  et  al.,  1988).  If  we  combine  these  closure  phases  with 
the  Fourier  magnitude  estimates  for  each  of  the  ( N“p )  baselines,  we  now  have  a  set  of 
atmospheric-invariant  observables  with  which  we  can  attempt  image  reconstruction. 
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The  Atmosphere  Problem 


Vjk  =  Vjk^i 


T 


Fourier  triple  products  ( bispectra )  cancel  atmosphere 

§123  =  V12V23V31 

ZV[2  =  012  +JP1  —  #2 
ZV23  =  023+1^  ->^3 
+  ^V3i  =  031  +jfi  ~ 


^§123  -  $12  +  $23  +  $31 

Figure  1.6:  Eliminating  the  effect  of  atmospheric  turbulence  with  phase  closure 
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Chapter  2 


Leveraging  compressed  sensing 
techniques  in  optical  interferometry 

2.1  Citation  to  Previously-Published  Work 

This  chapter  contains  text  and  figures  published  previously  in  the  following  paper: 

Kurien,  B.  G.,  Rachlin,  Y.,  Shah,  V.  N.,  Ashcom,  J.  B.  and  Tarokh,  V.  (2014).  Compressed 
sensing  techniques  for  image  reconstruction  in  optical  interferometry.  In  Imaging  and 
Applied  Optics  2014,  Optical  Society  of  America,  p.  SM2F.3. 

2.2  Chapter  Overview 

Optical  interferometers  image  a  scene  by  sampling  the  spatial  frequencies  that  comprise 
its  2D  Fourier  Transform.  The  sample  set  is  typically  much  smaller  than  the  number  of 
resolvable  elements  in  the  instrument's  field-of-view  (FOV).  In  the  related  field  of  radio 
interferometry,  sparse  recovery  (SR)  techniques  based  on  recent  seminal  work  (see,  e.g., 
Wiaux  et  al.  (2009))  have  proven  successful  in  regularizing  the  ill-posed  problem  arising 
from  this  undersampling.  In  contrast,  atmospheric  turbulence  precludes  the  availability  of 
direct  Fourier  phase  information  in  optical  interferometry,  and  hence  CS  techniques  do  not 
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directly  apply-  We  have  developed  and  validated  a  robust  algorithmic  interface  between  the 
Fourier  magnitude  and  bispectrum  observables  available  in  the  optical  interferometry  and 
SR  regularizations  known  as  Total  Variation  Minimization  for  which  fast  solvers  exist. 


2.3  Problem  Statement  and  Approach 

There  are  two  principal  challenges  that  one  encounters  when  trying  to  reconstruct  im¬ 
ages  from  optical  interferometric  measurements.  Each  aperture  pair  (or  baseline)  in  an 
interferometer  samples  the  Fourier  transform  of  the  image  in  its  field  of  view  at  a  spa¬ 
tial  frequency  j,  where  b  is  the  length  of  the  baseline,  and  A  is  the  wavelength  of  light 
collected.  Interferometers  rarely  measure  enough  baselines  to  fully  sample  the  Fourier 
transform.  Flence  reconstruction  of  the  image  from  such  undersampled  measurement  set  is 
an  ill-posed  (or  underdetermined)  problem;  we  must  recover  intensities  for  NreSEi  resolution 
elements  in  the  scene  under  observation  from  m  <C  Nn,SEl  spatial  frequency  measurements 
corresponding  to  the  available  baseline  pairs  in  the  telescope  array.  To  make  the  problem 
well-posed,  additional  constraints  must  be  enforced.  Traditionally  regularization  approaches 
based  on  the  so-called  CLEAN  algorithm  due  to  Hogbom  (1974),  which  implicitly  fits  the 
measurements  to  an  image  model  consisting  of  a  sparse  collection  of  point-like  sources, 
have  been  used  to  solve  this  problem. 

It  is  well-known  in  imaging  that  successful  regularization  schemes  apply  constraints 
which  accurately  capture  the  properties  of  the  image  we  wish  to  reconstruct.  Likewise  the 
success  of  the  CLEAN  algorithm  has  been  largely  limited  to  astronomical  scenes  which 
closely  match  its  point-source-collection  assumption.  The  astronomical  community  has 
hence  been  led  to  consider  more  generally-applicable  prior  models.  A  property  shared  by 
an  overwhelmingly-large  fraction  of  natural  images  is  compressibility.  Compressible  images 
are  those  which  can  be  well-approximated  by  a  sparse  representation  in  some  domain.  The 
smallness  of  an  image's  Ll-norm  1,  as  applied  to  either  directly  to  the  image  pixels,  their 

1The  Ll-norm  of  a  vectorized  image  x,  denoted  || x||  j,  is  the  sum  of  the  absolute  values  of  the  entries  in  x 
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gradient,  or  their  wavelet  coefficients,  has  proven  to  be  a  remarkably  effective  proxy  for  the 
image's  compressibility.  This  observation  and  techniques  that  exploit  it  are  at  the  core  of 
the  field  of  sparse  recovery  (SR)  as  well  as  the  tightly-linked  field  of  compressed  sensing  (CS) 
(Candes  et  al.,  2006)  (Donoho,  2006).  Though  a  wide  variety  of  sp arse-recovery  techniques 
have  arisen  in  the  past  few  decades,  the  vast  majority  of  them  regularize  problems  with  a 
linear  observation  model: 


y  =  Fx  +  n  (2.1) 

where  y  is  the  measurement  vector,  x  is  a  vector  representing  the  unknown  signal  or 
image,  F  is  the  measurement  matrix,  and  n  represents  additive  measurement  noise. 

While  sparse-recovery  techniques  have  been  applied  successfully  to  reconstruction 
in  radio  interferometry  (Wiaux  et  al.,  2009),  optical  interferometry,  on  the  other  hand, 
poses  additional  challenges  beyond  the  Fourier-undersampling  problem  discussed  above. 
Namely,  atmospheric  turbulence  necessitates  a  non-linear  formulation  of  the  reconstruction 
problem.  Recall  from  Chapter  1  that  there  are  ( N“p )  unknown  Fourier  phases  as  well  as 
Nnp  —  1  unknown  piston  differences,  and  hence  inference  of  the  Fourier  phases  from  the 
(N2P)  phase  measurements  available  in  a  non-redundant  array  is  inherently  ill-posed.  The 
traditional  means  of  mitigating  this  issue  has  been  to  again  impose  prior  constraints  on  the 
reconstructed  image  (e.g.  constraints  on  the  scene's  compactness,  sparsity,  or  smoothness). 
In  particular,  a  myriad  of  so-called  self-calibration  algorithms  have  been  developed  (see,  e.g., 
Pearson  and  Readhead  (1984))  which  alternate  between  estimation  of  the  Fourier  phases 
and  estimation  of  the  image  coefficients  subject  to  the  imposed  constraints.  Such  algorithms 
have  been  shown  to  be  successful  in  high-flux  scenarios  in  which  the  interferometric 
measurements  have  a  high  signal-to-noise  ratio  (SNR).  In  low-flux  scenarios  on  the  other 
hand,  one  must  integrate  observables  over  a  long  time  period  in  order  to  build  sufficient 
SNR  for  reliable  image  reconstruction.  Since  integration  beyond  the  coherence  time  of  the 
atmosphere  would  result  in  fringe  blurring,  one  is  then  led  to  the  formation  and  subsequent 
integration  of  atmosphere-invariant  observables  from  each  atmosphere-coherence  time  (or 
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frame). 

In  Section  1.6  we  introduced  the  bispectrum  observable  as  the  classical  basis  for 
atmosphere-invariant  inference  in  optical  interferometry.  Formed  as  the  product  of  three 
fringe  phasors  associated  with  the  sides  of  a  baseline  triangle  (e.g.  bi2,  b23,  and  b3i), 
the  bispectrum  is  a  non-linear  function  of  the  complex  visibilities  of  the  image.  Recall 
that  the  atmosphere-induced  terms  cancel  in  these  products  and  hence,  like  the  Fourier 
magnitudes,  these  so-called  bispectra  are  atmosphere-invariant  observables.  However,  for  a 
non-redundant  array  with  (^)  distinct  baselines,  recovery  of  the  Fourier  phases  from  the 
bispectra  phases  (i.e.  the  closure  phases )  remains  ill-posed  since  there  are  only  ( N2  1 )  indepen¬ 
dent  closure  phases  (Readhead  et  al,  1988).  Successful  bispectra-based  image  reconstruction 
remains  feasible  in  spite  of  this  ill-posedness  (see  e.g.  Thiebaut  (2013),  Besnerais  et  al.  (2008)), 
but  again  prior  constraints  (e.g.  on  the  image  support)  must  be  enforced  to  regularize  the 
reconstruction. 

An  important  property  of  the  bispectrum  which  we  will  exploit  in  our  reconstruction 
method  is  its  invariance  to  scene  translation.  To  illustrate  this  property,  let  us  denote  the 
spatial  frequency  vectors  of  two  sides  of  a  bispectrum  triangle  as  u  and  v,  respectively. 
Then  the  spatial  frequency  of  the  remaining  side  will  be  w  :=  —  (u  +  v).  If  y(p)  is  the 
corresponding  (normalized)  bispectrum  of  a  scene  at  its  original  position  p,  then  by  the 
shift  property  of  the  Fourier  transform,  the  bispectrum  of  the  translated  image  becomes: 

£(p_l_<5p)  =  ej{sn+u-Sp)ej{ev+vSp)ej(e_(u+v)-(u+v)-Sp)  (2.2) 

=  ei(e  u+0v+0_(u+v))  =  £(p)  (2.3) 

where  we  have  dropped  the  aperture  subscripts  for  purposes  of  generality. 

Since  the  bispectra  are  non-linear  functions  of  their  underlying  image,  reconstruction 
from  these  observables  does  not  map  directly  onto  the  linear  SR  framework  in  Equation  (2.1). 
To  address  this  challenge,  we  developed  an  interface  which  first  recovers  Fourier  component 
estimates  from  bispectra  and  Fourier  magnitudes,  using  a  nonlinear  least  squares  solver. 
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An  Algorithmic  Interface 


Figure  2.1:  Overview  of  Proposed  Two-Stage  Approach 


With  the  Fourier  component  estimates  at  hand,  the  image  estimate  can  then  be  recovered 
using  standard  SR  techniques.  The  approach  is  diagrammed  in  Figure  2.1  below. 


2.4  Description  of  Stage  1  of  Algorithm 

The  goal  of  Stage  1  of  our  algorithm  is  to  use  the  measured  bispectra  and  Fourier  magnitudes 
to  recover  estimates  of  the  true  Fourier  components  of  the  image.  Recall  from  Section  1.3 
that  the  fringes  provide  measurements  of  the  Fourier  magnitudes  directly.  Hence  it  remains 
to  recover  the  Fourier  phase  information.  Since  there  are  ( )  unknown  Fourier  phases  in 
a  non-redundant  array  but  only  (Na,f  ')  independent  closure  phases,  this  sub-problem  is 
itself  ill-posed.  Namely,  we  cannot  hope  to  find  a  unique  solution  for  the  Fourier  phases 
simply  by  minimizing  the  residuals  of  the  bispectra  data.  Instead  our  algorithm  employs  a 
gradient-based  search  to  find  the  phase  set  that  minimizes  a  regularized  objective  consisting 
of  a  data  term  and  a  prior  term.  The  data  term  is  the  sum  of  the  squared  bispectrum 
residuals.  The  prior  term  is  a  metric  that  favors  those  Fourier  phase  vectors  corresponding 
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to  a  spatially-compact  intensity  profile  in  the  image  domain.  This  kind  of  joint  metric  has 
been  suggested  for  optical  interferometric  applications  before  in  the  work  of  Thiebaut  (2013). 
However,  whereas  the  method  in  Thiebaut  (2013)  searches  for  a  metric-minimizing  image, 
our  method  first  searches  for  the  metric-minimizing  baseline  phase  set,  which  is  inherently 
a  smaller  set  given  the  undersampled  nature  of  the  problem. 


9  =  argminfdflte(0,g,a)  +  }tfprior{d,a)  (2.4) 

e 

The  data  (or  residual)  and  prior  terms  are  defined  as  follows. 


f  g,a) 


Nb 

L 


k=l 


1 

-  ( gk  ~  aklak2ak3exPij{6kl  +  ®k2  +  fys))) 

lVk 


(2.5) 


where  gk  is  the  k- th  bispectrum  measurement,  {6ki}  are  the  unknown  Fourier  phases  in 
the  k- th  bispectrum,  the  {flja}  are  estimates  for  the  magnitudes  of  the  Fourier  coefficients  in 
the  k- th  bispectrum,  and  ay  is  a  weighting  factor  proportional  to  the  estimated  variance  of 
the  fc-th  bispectrum. 


f prior (0,  a)  =  ||Ke[F*ya/0]  o  (1  -  h.gaussian)\\2  (2.6) 

where  F  is  the  partial  DFT  matrix  whose  rows  are  2D-sinusoidal  basis  functions  sampled 
by  the  array,  ya/e  is  the  vector  of  estimated  Fourier  coefficients  with  magnitudes  given  by  a 
and  phases  given  by  9,  h gaussian  denotes  a  centered  Gaussian  window  with  a  peak  value  of 
unity,  1  denotes  the  vector  whose  entries  are  all  1,  and  the  operation  o  denotes  point- wise 
multiplication. 

An  intuitive  interpretation  of  the  regularization  above  can  be  obtained  by  decomposing 
the  computation  into  steps.  We  seek  to  fit  the  bispectra  with  a  set  of  Fourier  phases  9  which 
together  with  the  estimated  Fourier  magnitudes  produce  an  image  with  compact  energy. 
The  data  term  in  the  metric  f ^ata  is  a  sum  of  the  squared  residuals  of  measured  bispectra 
with  respect  to  the  bispectra  associated  with  the  phase  vector  9.  The  steps  for  computation 
of  the  prior  term  are  as  follows.  First  we  associate  the  phases  in  9  with  their  corresponding 
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Figure  2.2:  Inverted  Gaussian  penalty  function  (l  —  h  gaussian)l  the  dark  region  represents  an  area  of  near-zero 
penalty 

magnitude  estimates  in  a  to  form  estimates  for  the  complex  Fourier  components  y3/e.  We 
then  apply  the  pseudo-inverse  of  the  partial  DFT  matrix  F  to  the  vector  ya  g  and  retain  the 
real  part  2 ,  thereby  obtaining  an  image-space  representation  of  ya/g.  Finally,  we  perform 
point-wise  multiplication  of  this  image  with  a  function  of  the  form  shown  in  Figure  and 
compute  the  norm  of  the  result.  This  operation  computes  a  compactness  metric  which 
penalizes  y3/g  in  proportion  to  its  spatial  energy  spread.  Note  that  the  actual  position  of  the 
Gaussian  taper  is  immaterial;  as  established  above,  the  bispectrum  is  a  translation-invariant 
quantity  and  hence  the  joint  metric  will  be  invariant  to  shifts  of  the  same  object. 

The  joint  data-prior  metric  in  Equation  (2.4)  can  be  minimized  using  a  suitable  non¬ 
linear  least-squares  solver.  For  our  processing,  we  have  elected  to  use  the  MATLAB® 
implementation  (MATLAB,  2012)  of  the  trust-region  reflective  non-linear  least-squares 
solver. 


2.5  Description  of  Stage  2  of  Algorithm 

Stage  2  takes  the  Fourier  estimates  obtained  in  Stage  1  and  attempts  to  recover  an  image 
using  fast  sparse-recovery  (SR)  techniques.  This  amounts  to  solving  the  linear,  under¬ 
determined  inference  problem: 

2Since  the  rows  of  F  are  orthogonal,  the  pseudo-inverse  is  obtained  by  a  Hermitian  transpose  operation  F* 
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y  =  Fx  +  n 


(2.7) 


where  y  is  the  vector  of  Fourier  estimates  from  Stage  1,  F  is  the  partial-DFT  matrix  whose 
rows  are  the  complex-sinusoidal  basis  functions  corresponding  to  the  spatial  frequencies 
sampled  by  the  array3,  and  n  is  the  (unknown)  residual  error  in  these  estimates.  Since 
in  practice  there  are  fewer  available  Fourier  measurements  than  the  desired  number  of 
resolution  elements  in  the  image,  this  inference  problem  is  ill-posed  from  a  Nyquist-sampling 
perspective.  Hence  once  again  we  must  regularize  the  problem  by  imposing  constraints  on 
the  image.  We  selected  two  different  SR  regularization  strategies  due  to  their  widespread 
success  in  solving  similar  inverse  problems  involving  Fourier  undersampling.  The  first 
of  these  approaches.  Basis  Pursuit  Denoising  (BPDN),  seeks  to  find  the  image  of  smallest 
Ll-norm  which  also  agrees  with  the  estimated  Fourier  component  estimates  to  within  a 
certain  tolerance  e.  We  performed  this  Ll-minimization  in  the  wavelet  domain  as  opposed 
to  the  pixel  domain,  as  it  is  well-known  that  natural  images  tend  to  be  more  compressible 
in  the  former  domain  than  in  the  latter.  The  other  regularization  we  tested.  Total  Variation 
Minimization  (TV-Min),  is  very  similar,  but  seeks  instead  to  minimize  another  quantity 
known  to  be  compressible  in  natural  images:  the  Ll-norm  of  the  gradient  of  the  image.  The 
two  regularization  techniques  are  specified  mathematically  below: 


(Basis  Pursuit  Denoising):  x  =  argmin  | j a  1 1  ,  subject  to: 

a 


FY  -  y 


<  e 

2 


(2.8) 


(Total  Variation  Minimization):  x  =  argmin  || a ||Ty  subject  to:  |Fa  —  y  ||2  <  e  (2-9) 

a 


where  ||a||rv-  is  the  Ll-norm  of  the  2D  pixel  gradient  of  the  image,  i.e.  \\a\  TV 


E/,/HV4',/]||1. 


3Note  that  F  need  not  be  the  same  as  F,  as  we  are  free  to  assume  an  independent  pixel  resolution  in  each 
case.  This  resolution  then  specifies  the  granularity  of  the  sampling  of  the  sinusoids  forming  the  rows  of  the 
matrix. 
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The  software  package  NESTA  (Becker  et  al.,  2011)  was  used  to  perform  the  Stage  2 


optimizations  above. 


2.6  Algorithm  Performance 

2.6.1  Laboratory  Validation 

To  test  the  effectiveness  of  our  algorithm  when  applied  to  real  interferometric  data,  we 
collected  fringe  data  with  MIT  Lincoln  Laboratory's  (MIT/LL)  Lizeau  Interferometer.  The 
interferometer  is  a  reconfigurable  fiber-coupled  system  that  allows  all  (22°)  baseline  conjugate 
pairs  to  be  measured  simultaneously  by  projecting  light  from  all  apertures  onto  a  Lizeau 
beam-combiner.  Lor  our  experiments,  we  used  a  20-aperture  compact,  non-redundant 
pattern  of  the  Golay  type  (Golay,  1970).  The  pattern  is  shown  in  Ligure  2.3.  The  ratio 
r  of  the  maximum  baseline  to  the  minimum  baseline  was  17  in  both  x  and  y  spatial 
coordinates.  Hence  the  number  of  resolution  elements  required  to  uniquely  specify  our 
diffraction-limited  scene  was  (2 r  )2,  yielding  an  undersampling  ratio  of  p  =  ~  0.33. 

The  interferometer  was  illuminated  with  the  far-field  projection  of  a  range  of  targets:  chrome- 
on-glass  transparency  masks  were  photolithographically-prepared  and  placed  at  the  focus 
of  an  j2  three-mirror  off-axis  telescope.  When  the  target  is  illuminated  with  the  white  light, 
the  far-field  projection  of  the  target  is  produced  at  the  telescope  aperture,  where  it  can  be 
sampled.  As  a  reference  for  reconstruction,  we  consider  a  reduction  of  the  actual  target  to 
the  fundamental  diffraction-limited  resolution  of  the  instrument  as  shown  in  Ligure  2.4. 

5000-frame  reconstruction  results  are  shown  in  Ligure  2.5  for  the  following  photoelectron 
(pe)  levels:  (from  left  to  right)  320<?3  pe/ aperture /frame,  54<?3  pe/ aperture /frame,  and  1 5e3 
pe/ aperture/ frame. 

2.6.2  Simulation 

In  addition  to  the  Laboratory  validation,  we  also  conducted  5000-frame  simulations  with 
the  same  aperture  pattern  to  compare  algorithm  performance  in  the  presence  of  unknown 
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Aperture  Pattern  Frequency  Sampling 
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Figure  2.3:  N on-redundant  Golay  20-aperture  pattern  (left)  and  corresponding  UV-sampling  (right) 


>< 


Figure  2.4:  Target  chrome  mask  (left)  and  corresponding  truth  image  at  diffraction-limited  resolution  (right) 
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Figure  2.5:  Image  Reconstruction  Results  from  Laboratory  Validation 
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Figure  2.6:  Image  Reconstruction  Results  from  Simulations 


atmospheric  distortion  against  that  if  this  distortion  were  known.  Note  that  in  the  latter 
atmosphere-oracle  case,  bispectrum  formation  is  not  required  and  direct  SR  reconstruction 
from  complex  amplitudes  is  feasible.  Hence  this  comparison  allows  a  quantitative  measure  of 
the  reconstruction  penalty  suffered  from  bispectrum  formation.  Our  comparative  simulation 
scheme  is  depicted  in  the  block  diagram  in  Figure  2.6. 

We  used  two  Image  Quality  Metrics  to  assess  performance:  the  standard  metric  Normal¬ 
ized  Mean-Squared  Error  (NMSE)  and  the  Structural  Similarity  (SSIM)  metric  (Wang  et  al., 
2004).  The  NMSE  is  defined  as: 


NMSE  =  20  log  ^ ^  (2.10) 

x0 

Neither  metric  captured  the  visual  similarity  of  the  raw  reconstructions  of  our  two-stage 
algorithm  relative  to  ground  truth.  To  address  this,  we  uniformly  hard-thresholded  each 
of  the  raw  reconstruction  pixels  before  evaluating  the  metrics.  Figure  2.7  shows  our  two- 
stage  reconstructions  before  thresholding  (left),  the  same  reconstructions  after  thresholding 
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Figure  2.7:  Image  Reconstruction  Results  from  Simulations 
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Figure  2.8:  Algorithm  Performance  in  Simulation 


(middle),  and  the  atmosphere-oracle  reconstructions  (right).  Results  from  both  a  low-light- 
level  scenario  (top  row)  and  a  high-light-level  scenario  (bottom  row)  are  shown.  Figure  2.8 
uses  NMSE  and  SSIM  image  quality  metrics  to  quantify  reconstruction  performance  for  our 
two-stage-algorithm  reconstructions  after  thresholding  (red),  and  the  atmosphere-oracle 
reconstruction  (black).  We  see  that  the  algorithm  achieves  performance  close  to  that  of  the 
atmosphere  oracle  for  a  wide  range  of  light  levels. 
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Chapter  3 


Pattern  design  criteria  for  uniqueness 
in  phase  recovery 

3.1  Citation  to  Work  under  Review 

This  chapter  contains  text  and  figures  submitted  for  publication  to  the  Monthly  Notices  of 
the  Royal  Astronomical  Society  (MNRAS). 

Kurien,  Bv  Tarokh,  V.,  Ashcom,  J.,  Rachlin,  Y.  and  Shah,  V.  (2016).  Resolving  phase 
ambiguities  in  the  calibration  of  redundant  interferometric  arrays:  implications  for  array 
design.  Monthly  Notices  of  the  Royal  Astronomical  Society,  (submitted,  March  4,  2016) 


3.2  Chapter  Overview 

We  provide  new  results  enabling  robust  interferometric  image  reconstruction  in  the  presence 
of  unknown  aperture  piston  variation  via  the  technique  of  Redundant  Spacing  Calibration 
(RSC).  The  RSC  technique  uses  redundant  measurements  of  the  same  interferometric 
baseline  with  different  pairs  of  apertures  to  reveal  the  piston  variation  among  these  pairs. 
In  both  optical  and  radio  interferometry,  the  presence  of  phase-wrapping  ambiguities  in 
the  measurements  is  a  fundamental  issue  that  needs  to  be  addressed  for  reliable  image 
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reconstruction.  In  this  paper,  we  show  that  these  ambiguities  affect  recently-developed 
RSC  phasor-based  reconstruction  approaches  operating  on  the  complex  visibilities,  as  well 
as  traditional  phase-based  approaches  operating  on  their  logarithm.  We  also  derive  new 
sufficient  conditions  for  an  interferometric  array  to  be  immune  to  these  ambiguities  in 
the  sense  that  their  effect  can  be  rendered  benign  in  image  reconstruction.  This  property, 
which  we  call  wrap-invariance,  has  implications  for  the  reliability  of  imaging  via  classical 
three-baseline  phase  closures  as  well  as  generalized  closures.  We  show  that  wrap-invariance 
is  conferred  upon  arrays  whose  interferometric  graph  satisfies  a  certain  cycle-free  condition. 
For  cases  in  which  this  condition  is  not  satisfied,  a  simple  algorithm  is  provided  for 
identifying  those  graph  cycles  which  prevent  its  satisfaction.  We  apply  this  algorithm  to 
diagnose  and  correct  a  member  of  a  pattern  family  popular  in  the  literature. 

Before  we  discuss  phase-wrapping  ambiguities,  we  first  review  a  few  key  preliminary 
mathematical  notions  from  lattice  theory  in  the  next  Section. 

3.3  Preliminaries 

3.3.1  Lattices 

A  lattice  is  a  mathematical  object  describing  a  repeating  pattern  of  discrete  points  in  space. 
It  arises  in  a  variety  of  scientific  and  engineering  contexts  including  crystallography  and 
communication  theory.  We  will  define  a  lattice  A  as  a  linear  combination  of  vectors  in  which 
the  coefficients  are  integers,  i.e. 

flfV/  |  Vaf  e  Z  j  (3.1) 

If  the  generating  vectors  v,  are  not  linearly-independent,  one  might  ask  for  another, 
more  compact  description  of  A.  Consider  one  such  candidate  set  of  linearly-independent 
vectors  b;.  If  all  integer  linear  combinations  of  b;  lie  in  A,  and  every  lattice  point  in  A  can 
be  represented  as  an  integer  linear  combination  of  the  by,  then  we  say  that  the  vectors  v; 
form  a  basis  for  the  lattice  A. 
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Many  lattice  algorithms,  including  algorithms  for  the  Closest  Vector  Problem  introduced 
in  the  next  section,  require  a  reduced  lattice  basis.  A  reduced  basis  is  one  consisting  of  vectors 
which  are  short  and  nearly-orthogonal.  A  well-known  procedure  to  create  a  basis  for  a 
given  lattice  with  short,  nearly-orthogonal  vectors  is  known  as  the  LLL  algorithm  due  to 
Lenstra  et  al.  (1982). 

An  LLL-based  routine  for  forming  a  reduced  lattice  basis  from  a  set  of  generating  vectors 
which  are  not  necessarily  linearly-independent  is  given  in  Cohen  (1993). 


3.3.2  The  Closest- Vector-Problem 

A  well-known  problem  in  lattice  theory  is  the  Closest-Vector-Problem,  which  can  be  described 
as  follows: 

Problem  3.3.1.  Given  a  basis  {bi,b2,  ...,b n}  for  a  lattice  in  Rm  and  a  vector  w  G  Rm,find  the 
point  in  the  lattice  closest  in  Euclidean  distance  from  w. 

Several  algorithms  exist  for  solving  this  problem.  A  popular  class  of  algorithms,  known 
as  the  Sphere-Decoding  algorithms,  are  efficient  searches  for  the  closest  lattice  point  within 
a  hypersphere  of  a  certain  radius  centered  on  the  input  vector  (see  e.g.  Agrell  et  al.  (2002)). 
For  the  simulations  in  this  chapter,  we  instead  use  the  lower-complexity  Babai  Nearest 
Plane  (Babai-NP)  algorithm  (Babai,  1986).  For  lattice  bases  which  are  nearly  orthogonal,  this 
algorithm  offers  reliable,  albeit  not  guaranteed,  performance  in  practice.  Pseudo-code  for 
one  implementation  of  this  algorithm  due  to  Galbraith  (2012)  is  given  in  listing  Algorithm  1. 


Algorithm  1  Babai  Nearest-Plane  Algorithm 


Input:  Basis  for  lattice  A  (i.e.  {bi, . . .  ,b„}),  and  w  G  1R”' 

Output:  Element  v*  nearest  to  w  in  A 

Compute  orthogonal  basis  {bp  . . . ,b* }  using  Gram-Schmidt  procedure 
for  i  =  n  downto  1  do 


Set  /•  - 
Set  y i  =  [il}bl 

Set  w/_i  =  W;  -  (U  -  [li])b*  -  [li]bi 


end  for 

return  v  =  yi  +  . . .  y„ 
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A  conceptual  interpretation  of  this  recursive  algorithm  is  as  follows.  At  the  root  level, 
the  algorithm  finds  the  vector  in  A  closest  to  w;  in  other  words  it  solves  the  CVP  problem 
(An,  w).  First  consider  the  subspace  U  spanned  by  the  first  n  —  1  basis  vectors,  i.e.  U  = 
span{ bi, . . .  The  algorithm  begins  by  computing  the  translate  y  G  A  of  U  closest 

to  w,  i.e.  the  nearest  plane  to  w.  The  algorithm  then  implicitly  computes  the  projection 
w'  of  w  onto  this  translate  U  +  y.  It  then  translates  this  projection  back  to  the  origin 
(i.e.  w„_i  :=  w'  —  y),  and  solves  the  smaller  problem  within  U,  i.e.  the  CVP  problem 
(A„_i,  wn_i),  where  A„_i  is  the  lattice  formed  by  the  first  n  —  1  basis  vectors.  The  process 
continues  recursively  until  the  final  CVP  problem  with  the  lattice  spanned  by  bi  is  solved. 
The  final  output  is  then  the  sum  of  the  translates  from  each  level  of  the  recursion. 

3.4  Problem  Statement  and  Related  Work 

We  have  seen  that  unknown  optical  path  differences  (OPD)  amongst  the  apertures  in  an 
array,  arising  from  atmospheric  turbulence  as  well  as  non-idealities  in  the  interferometric 
system,  present  a  fundamental  challenge  in  interferometry.  Moreover,  we  have  seen  that 
while  the  bispectrum  observable  eliminates  the  OPD  in  interferometric  measurements, 
the  ill-posed  nature  of  Fourier  phase  recovery  arising  from  the  OPD  terms  persists.  In 
Chapter  2  we  developed  an  algorithmic  interface  which  regularizes  the  ill-posed  problem  of 
recovering  Fourier  phase  from  the  bispectrum  by  enforcing  prior  constraints  on  the  image. 
An  alternative  approach  to  prior-regularized  reconstruction  is  to  use  baseline  redundancy 
to  explicitly  solve  for  OPD  variation;  an  array  with  baseline  redundancy  contains  repeated 
instances  of  the  same  vector  baseline  involving  distinct  aperture  pairs.  Since  Fourier  phases 
can  be  assumed  to  be  equal  for  all  repeated  vector  baselines,  an  observed  difference  amongst 
their  corresponding  measurements  exposes  the  contribution  of  the  OPD.  This  idea  of  using 
redundant  arrays  to  calibrate  out  OPD  variation,  known  as  redundant  spacing  calibration 
(RSC),  was  developed  in  works  such  as  those  by  Arnot  et  al.  (1985)  and  Greenaway  (1990). 
In  recent  years,  innovation  in  optical  technology  has  engendered  a  revival  of  interest  in  the 
RSC  technique.  The  simultaneous  (or  Fizeau- style)  measurement  of  fringes  on  a  common 
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Figure  3.1:  Fraction  of  redundant  baselines  required  for  critical  redundancy  vs.  aperture  count 

focal  plane  has  long  been  a  popular  method  of  acquiring  many  baseline  measurements 
in  an  economical  manner.  However,  the  Fizeau  method  had  been  incompatible  with  RSC 
techniques  since  the  fringes  formed  by  each  set  redundant  baselines  would  alias  on  the 
focal  plane.  An  elegant  solution  to  this  problem  was  proposed  by  Perrin  et  al.  (2006).  This 
work  developed  the  idea  of  segmenting  the  entrance  pupil  of  a  single  telescope  into  an 
RSC  arrangement  of  sub-pupils  from  which  the  light  was  then  coupled  via  single  mode 
fiber  to  a  non-redundant  exit  pupil,  thereby  permitting  unambiguous  and  simultaneous 
fringe  detection  for  an  RSC  array.  A  reconstruction  algorithm  for  this  architecture  was 
then  proposed  in  Lacour  et  al.  (2007).  Even  more  recently,  RSC  has  been  implemented  as 
the  calibration  scheme  of  choice  for  several  radio  interferometers:  the  Donald  C.  Backer 
Precision  Array  for  Probing  the  Epoch  of  Reionization  (PAPER)  in  South  Africa  (see  Ali 
et  al.  (2015)),  the  MIT  Epoch  of  Reionization  (MITeOR)  in  the  United  States  (see  Zheng  et  al. 
(2014)),  and  the  Ooty  Radio  Telescope  (ORT)  in  India  (see  Marthi  and  Chengalur  (2014)). 

As  will  be  shown  below,  N  —  3  independent  redundant  relations  are  required  for  unique 
determination  of  atmosphere  and  Fourier  phases  -  a  condition  we  will  refer  to  as  critical 
redundancy.  An  oft-cited  drawback  of  the  RSC  approach  is  that  it  reduces  the  number  of 
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unique  spatial  frequencies  measured  by  the  interferometer.  However,  as  Figure  3.1  illustrates, 
the  fraction  of  distinct  uv-samples  sacrificed  for  critical  redundancy  becomes  increasingly 
negligible  as  the  number  of  apertures  in  the  array  increases.  Nevertheless  the  RSC  technique 
presents  other  challenges  which  must  be  overcome  for  reliable  imaging.  Central  among 
these  challenges  is  the  problem  of  integer  phase  ambiguities  which  arise  from  the  fact  that 
the  interferometric  phase  is  only  known  modulo  2n.  Indeed  these  ambiguities  have  been 
shown  to  play  an  important  role  in  accurately  recovering  sensor  complex  gains  and  object 
complex  visibilities  while  imaging  with  real  interferometric  instruments  (see  e.g.  Liu  et  al. 
(2014),  Eastwood  et  al.  (2009)).  In  this  chapter,  we  describe  these  ambiguities  and  how  they 
can  be  mitigated  using  a  combination  of  lattice  theory  algorithms  and  careful  array  design. 
We  will  see  that  these  ambiguities  have  a  fundamental  presence;  namely,  they  exist  whether 
the  calibration  strategy  works  with  complex  visibilities  (which  we  call  the  Phasor  approach) 
directly,  or  their  respective  logarithms  (which  we  call  the  Phase  approach).  To  the  best  of 
our  knowledge,  the  results  in  this  chapter  are  the  first  to  provide  array  conditions  allowing 
unambiguous  interferometric  phase  determination  in  spite  of  wrap  ambiguities  in  the  Phase 
approach,  and  corresponding  false  minima  in  the  objective  of  the  Phasor  approach. 

To  motivate  the  analysis  in  the  chapter,  we  provide  an  example  illustrating  the  effect 
that  wrap  ambiguities  can  have  in  RSC-based  image  reconstruction.  Consider  the  pattern 
depicted  in  the  Figure  3.2.  This  pattern  belongs  to  one  of  the  more  popular  array  classes  in 
the  interferometry  literature:  the  so-called  Y-patterns  (see  e.g.  Arnot  et  al.  (1985),  Blanchard 
et  al.  (1996),  Labeyrie  et  al.  (2006),  Eastwood  et  al.  (2009),  Liu  et  al.  (2014)).  The  corresponding 
spatial,  or  UV,  sampling  is  provided  in  the  right  panel  of  the  Figure. 

To  demonstrate  the  potential  effect  of  wrap  ambiguities  on  reconstruction,  we  simulated 
both  Phase-  and  Phasor-based  reconstruction  from  noiseless  measurements  with  this  pattern. 
The  results  are  shown  in  Figure  3.3.  The  upper  left  panel  shows  the  true  image,  and  the 
upper  right  panel  shows  the  inverse  Fourier  transform  of  the  UV-sampled  visibility  function 
of  the  object  (i.e.  the  so-called  dirty,  or  interferometric  image).  The  lower  left  panel  shows  the 
reconstruction  result  with  an  implementation  of  the  Phase  method  similar  to  that  developed 
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Y-Type  Aperture  Pattern  UV-Sampling 


Figure  3.2:  Y-pattern  Array  Example 

by  Lannes  (2003),  and  the  lower  right  panel  shows  the  same  for  an  implementation  of  the 
Phasor  method  (Marthi  and  Chengalur,  2014)  (Lacour  el  al.,  2007).  Reconstruction  suffers 
from  phase  wrapping  error  in  the  former  case,  and  a  corresponding  false-minimum  trap 
in  the  Phasor  case.  In  the  course  of  this  chapter,  we  will  first  identify  this  ambiguity  from 
a  mathematical  perspective,  relate  it  to  a  particular  physical  structure  (i.e.  the  existence 
of  a  certain  type  of  loop  in  the  interferometric  graph),  and  provide  a  simple  algorithm  for 
identifying  such  structures  in  an  arbitrary  array  so  that  they  can  be  remedied. 

The  chapter  is  organized  as  follows.  In  Section  3.5,  we  review  previous  work  on  the 
integer  ambiguity  problem,  and  discuss  its  presence  in  the  Phase  approach.  We  provide 
new  mathematical  conditions  for  an  aperture  pattern  to  be  wrap-invariant ,  meaning  that 
the  effect  of  the  271-periodicity  of  the  measured  interferometric  phase  can  be  eliminated 
in  image  reconstruction.  These  results  are  founded  upon  techniques  from  lattice  theory, 
as  well  as  the  well-known  Smith  Normal  Form  (SNF)  of  an  integer  matrix.  We  show  the 
implications  of  these  results  on  imaging  with  three  types  of  interferometric  observables: 
the  baseline  phase  measurements,  their  traditional  closure  phases,  and  generalized  closure 
phases.  In  Section  3.6,  we  relate  these  mathematical  conditions  to  conditions  on  the  aperture 
pattern  itself.  Namely  we  show  that  wrap-invariance  is  conferred  upon  arrays  satisfying  a 
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Phase  Approach  (Noiseless)  Phasor  Approach  (Noiseless) 


Figure  3.3:  Reconstruction  Results  for  Y-Pattern  Example 
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certain  loop-free  condition.  As  an  illustrative  example,  we  diagnose  a  pattern  belonging 
to  the  popular  Y-pattern  class  and  remedy  it  to  be  loop-free.  In  Section  3.7,  we  show  that 
the  computationally-complex  SNF-based  approach  for  ambiguity  resolution  is  actually  not 
necessary  for  wrap-invariance  in  many  cases,  as  long  as  a  wrap-induced  image  shift  can  be 
tolerated  1 .  For  such  scenarios,  we  show  that  wrap-invariance  provides  a  certificate  for  the 
success  of  existing  Phase  and  Phasor-based  approaches  in  avoiding  wrap-induced  errors 
in  the  former,  and  false  global  minima  in  the  latter.  Finally  we  summarize  our  results  in 
Section  3.8. 


3.5  Phase  Wrapping  Ambiguities  in  RSC  Image  Reconstruction 

In  this  Section,  we  describe  a  Phase  Approach  algorithm  leveraging  the  CVP-approach  for 
phase-wrap  resolution,  which  to  the  best  of  our  knowledge  was  first  developed  in  Lannes 
and  Anterrieu  (1999).  We  begin  by  identifying  the  fundamental  phase  ambiguity,  and 
subsequently  assess  its  impact  on  the  phase  error  in  the  RSC  phase  solution.  In  the  process, 
we  develop  mathematical  conditions  for  wrap-invariance  which  will  form  the  basis  for  the 
notion  of  a  wrap-invariant  pattern  in  Section  3.7.  Finally,  we  assess  our  results  in  the  context 
of  similar  results  for  approaches  requiring  the  use  of  closure  phases  (Lannes,  2003). 

3.5.1  Identifying  the  Fundamental  Phase  Ambiguity 

The  traditional  approach  to  RSC  reconstruction  operates  on  the  measured  baseline  phases 
(see  e.g.  Arnot  et  al.  (1985),  Greenaway  (1994)).  To  illustrate  the  approach,  let  us  consider  an 
interferometer  which  operates  at  a  wavelength  A  with  two  apertures  (say,  i  and  j)  separated 
by  a  vector  baseline  distance  of  b;;.  In  the  absence  of  any  optical  path  difference,  the 
interference  pattern  formed  by  these  two  apertures  encodes  a  sample  of  the  object's  Fourier 
Transform  at  spatial  frequency  Let  the  true  Fourier  phase  (which  we  will  refer  to  as 

*As  will  be  shown,  this  image  shift  is  distinct  from  the  fundamental  degeneracy  of  object  position  in 
interferometric  measurements 
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object  phase),  measured  by  this  interference  pattern  be  denoted  as  Op.  The  measured  phase  is 
given  by: 


fiij  —  Op  +  ( pj  —  (pi  +  2ne  (3.2) 

where  (pj  —  (pi  is  the  optical  path  difference  between  apertures  j  and  i,  and  e  is  unknown 
phase  wrap  integer  arising  from  the  fact  that  interferometric  phase  measurements  are  only 
known  modulo  2n. 

Consider  an  interferometric  array  which  simultaneously  makes  many  such  measurements 
amongst  its  N  apertures.  Suppose  that  of  the  array's  ( (j )  baselines,  d  of  them  are  distinct. 
Further  suppose  we  have  a  solution  set  {(pi}  and  { Op }  for  these  equations.  Let  r,-  denote  the 
vector  position  of  the  z-th  aperture.  As  noted  by  several  authors  (see,  e.g.  Wieringa  (1992)), 
we  can  obtain  another  valid  solution  set  simply  by  replacing  each  (pj  with  (p{  =  (p,  +  <po  +  z  •  r ,, 
and  each  Op  with  0-’  =  Op  —  z  •  (tj  —  r ,-),  for  arbitrary  (po  and  z.  Since  the  free  vector  z  is  a 
two-parameter  vector  representing  the  inherently-ambiguous  position  of  the  image  within 
the  Field-of-View  and  the  free  parameter  (po  is  simply  a  scalar  piston  offset,  the  kernel  of  the 
RSC  system  is  three-dimensional.  This  is  the  tilt-position  degeneracy  which  is  fundamental  in 
interferometric  reconstruction.  The  end  result  is  that  the  RSC  system  contains  d  unknown 
distinct  object  phases  and  N  unknown  aperture  pistons,  and  is  rank-deficient  by  at  least 
3.  This  implies  that  there  are  at  most  (d  +  N  —  3)  linearly-independent  equations  in  the 
RSC  system,  and  hence  at  most  N  —  3  redundant  relations  can  be  linearly-independent. 
Remaining  relations  can  be  expressed  as  linear  combinations  of  the  elements  in  this  set. 
A  commonly-used  simple  example  of  such  dependencies  (Greenaway,  1994)  is  shown  in 
Figure  3.4.  The  phase  measurement  ^34  associated  with  the  baseline  bo, 4  can  be  expressed  as 
a  linear  combination  of  those  associated  with  the  other  baselines,  i.e.  4834  =  /I12  —  /l  1 3  +  p>24- 

We  will  assume  for  the  remainder  of  the  chapter  that  our  array  contains  N  —  3  indepen¬ 
dent  relations.  Under  this  assumption,  we  could  in  principle  solve  for  a  particular  solution 
of  this  system  by  arbitrarily  setting  two  object  phases  (whose  spatial  frequencies  are  not 
co-linear)  and  one  piston  phase.  This  particular  solution  would  then  differ  from  the  true 
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Figure  3.4:  Example:  Dependent  Parallelogram  redundancy. 


Figure  3.5:  Example:  5-aperture  RSC  pattern.  The  six  distinct  baselines  are  shown. 

solution  by  a  phase  ramp  in  the  Fourier  domain,  corresponding  to  an  image  shift  in  the 
spatial  domain.  In  this  section  we  will  instead  construct  a  different  particular  solution  which 
is  immune  to  the  effects  of  phase-wrapping  by  design. 

As  an  example,  consider  the  simple  array  in  Figure  3.5.  There  are  (®)  =  10  baselines,  of 
which  4  are  redundant.  A  critical  array  of  5  apertures  would  have  2  redundancies.  Therefore 
this  array  possesses  more  redundancies  than  necessary  (call  it  strongly  redundant),  and  we 
anticipate  that  the  resulting  system  will  be  overdeter  mined. 

The  measurement  equations  associated  with  this  array  can  be  written  in  matrix  form: 
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(3.3) 


Denoting  the  matrix  above  as  M,  we  can  write  such  systems  in  compact  form  as: 


M 


e 

(p 


=  /3  +  27re 


(3.4) 


The  example  above  illustrates  that  the  general  phase  measurement  matrix  will  have  two 
sets  of  columns:  one  corresponding  to  the  object  phases,  and  one  corresponding  to  the  path 
differences.  We  can  divide  the  range  of  the  matrix  into  two  respective  subspaces  as  follows: 


Definition  3.5.1.  Two  fundamental  subspaces  (Lannes  and  Anterrieu,  1999):  The  spectral 
phase  space  K  is  the  span  of  the  d  columns  associated  with  the  Fourier  phases,  or  equivalently 
the  range  of  the  sub-matrix  M(j  formed  by  these  columns.  The  piston,  or  aberration,  phase 
space  L  is  the  span  of  the  N„p  columns  associated  with  the  piston  phases,  or  equivalently 
the  range  of  the  sub-matrix  formed  by  these  columns. 

If  we  let  n  =  (^)  be  the  number  of  baselines  in  the  array,  the  phase  measurement  matrix 
M  will  be  of  size  n  x  (d  +  N).  For  a  strongly-redundant  array  like  the  one  in  the  example 
above,  the  column-space  K  +  L  of  the  matrix  will  clearly  not  span  R”.  Therefore  the  wrapped 
measurement  vector  /!  will  not  in  general  fall  in  the  the  subspace  K  +  L  (and  potentially 


45 


can  be  quite  far  from  it).  In  the  absence  of  measurement  noise,  we  can  unwrap  these 
measurements  by  identifying  those  integer  correction  vectors  e  for  which  fi*  =  /3  +  2rre 
lies  in  K  +  L.  In  the  presence  of  noise,  on  the  other  hand,  the  unwrapped  vector  will  not 
generally  lie  in  K  +  L  (but  for  low-to-moderate  noise  will  be  in  the  vicinity).  Thus  we  search 
for  vector(s)  fi*  which  are  as  close  to  K  +  L  as  possible  in  a  weighted  least-squares  sense 
(Lannes  and  Anterrieu,  1999),  i.e.  we  search  for  the  vector 


0RSC 

W 

( 

e 

\ 

tr  sc  = 

Qrsc 

=  argmine>04 

/T(e)  —  M 

<P 

) 

where  W  is  the  weighting  matrix.  If  we  let  E  denote  the  phase  measurement  covariance 
matrix  and  set  W  =  E  this  is  equivalent  to  searching  for  vectors  e  which  minimize  the 

projection  of  a  whitened  measurement  W2 fi*  =  W2  (/  +  2zre)  onto  the  space  (. K  +  L)w  := 

1  T 

ker( (WzM)  ),  where  ker(A)  denotes  the  kernel,  or  nullspace,  of  the  matix  A. 

Specifically  we  seek  to  minimize: 

/(e)  =  ||PwW2Q3  +  27re)||2  (3.6) 

where  Pw  is  a  matrix  representing  the  orthogonal  projection  from  IR"  onto  ( K  +  L)w. 
Letting  e'  =  — e,  we  can  rewrite  the  above  objective  function  as: 

2  2 

f(e')  =  ||PwWz(j3  — 27re')||  =  |  |PwWz/5  -  27rPwW2e'|  |  (3.7) 

Lannes  and  Anterrieu  (1999)  showed  that  this  optimization  problem  is  equivalent  to  the 
so-called  closest  vector  problem  (CVP)  in  the  theory  of  lattices.  We  will  define  a  lattice  L(Z") 
as  the  set  of  points  generated  by  integer  combinations  of  the  column  vectors  of  a  matrix 
L.  Letting  P  =  PwW?,  our  optimization  problem  then  amounts  to  the  following:  Find  the 
lattice  point  in  P(Z")  which  is  closest  to  P/3.  Indeed  this  CVP  formulation  applies  in  a  very 
general  sense  to  phase-unwrapping  problems  taken  collectively  as  discussed  in  Wubben 
et  al.  (2011). 

A  compact  representation  of  the  lattice  T  is  given  by: 
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(3.8) 


{m<n—(d+N— 3)  'l 

fl,v,  |  Vflf  ezj 

where  {v,}  are  linearly-independent  and  together  form  a  basis  of  the  lattice. 

Suppose  we  have  found  a  basis  for  the  lattice  P(Zn),  and  we  have  solved  the  Closest 
Vector  Problem  for  a  given  measurement  vector  /3  (c.f.  Section  3.3.2).  Let  b  be  the  output 
of  the  Babai  Nearest  Plane  Algorithm  -  i.e.  it  is  the  lattice  point  which  is  the  closest  to  jfb 
We  now  seek  to  solve  for  the  wrap  vector  corresponding  to  this  lattice  point,  i.e.  we  seek  a 
solution  to: 


b*  =  Pe  (3.9) 

Note  that  P  is  a  (weighted)  projection  matrix  and  thus  not  full-rank,  and  therefore  there 
will  be  infinitely-many  solutions  to  this  equation.  The  Closest- Vector-Problem  algorithm 
will  provide  one  particular  solution  ep.  The  complete  set  of  solutions  is  then  given  by: 

e  =  ep  +  eh  (3.10) 

where  ej,  is  any  integer  vector  in  the  kernel  of  P.  Suppose  we  choose  one  such  vector  e/., 
and  correct  our  phase  measurement  vector  accordingly.  The  corrected  phase  measurement 
vector  can  be  written  as: 

K  =  p  +  2n(ep  +  eh)  (3.11) 

Lemma  3.5.1.  e?,  G  K  +  L,  Ve;, 

Proof:  The  fact  that  e;,  G  ker( P^Wi)  implies  that  W z e /,  G  ker( Pw)-  This  in  turn  implies 
that  W^e/,  G  z’m(WzM)  and  hence  that  e/,  G  im( M)  since  is  invertible  by  construction. 
□ 

While  any  choice  of  e/,  G  K  +  L  will  admit  a  solution  to  Equation  (3.5),  let  us  consider 
the  optimal  one  e/,  o  which  minimizes  the  error  in  the  ultimate  phase  solution  Trsc-  Let 
£o  =  iS  +  27T(ep  +  e;,  o)  be  the  corresponding,  optimal  unwrapped  measurement  vector. 
From  Lemma  3.5.1,  we  see  that  the  unwrapped  vector  /I  differs  by  some  2ueJ1  in  K  +  L 
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from  j8q,  i.e. 


P  =  (i*+2nel  (3.12) 

The  impact  of  this  fundamental  ambiguity  is  the  main  subject  of  this  chapter.  We  have 
depicted  the  situation  in  Figure  3.6  to  provide  a  visual  summary  of  the  linear  algebra 
involved.  As  we  have  seen,  the  RSC  process  begins  with  the  interferometric  phase  mea¬ 
surement  /3,  which  due  to  wrapping  will  in  general  lie  far  from  the  range  K  -|-  L  of  the 
measurement  matrix.  By  solving  the  Closest- Vector-Problem  using  a  lattice  algorithm  such 
as  the  Babai  Nearest  Plane  algorithm,  it  is  possible  to  find  a  particular  correction  vector 
2nep  which  when  added  to  /3  minimizes  the  residual  in  Equation  (3.5),  i.e.  the  weighted 

/V  >(c 

distance  from  K  +  L.  In  the  noiseless  case,  this  unwrapped  vector  /3  will  lie  in  K  +  L  (i.e. 
zero  residual),  whereas  in  the  noisy  case,  it  will  in  general  remain  outside  of  K  +  L.2.  In 
either  case,  the  choice  of  the  residual-minimizing  vector  /T  is  not  unique.  To  see  this,  let 
the  smallest  possible  residual  norm  among  all  unwrapped  candidates  be  denoted  as  rmi„. 
The  set  of  unwrapped  measurement  vectors  rm;n  away  from  K  +  L  can  be  represented  as 
discrete  points  in  a  plane  parallel  to  K  +  L,  each  of  which  corresponds  to  a  distinct  choice  of 
ep.  Within  this  family,  consider  the  optimum  vector  fi(j  whose  corresponding  least-squares 
solution  r rsc  has  minimal  error.  Any  choice  of  an  unwrapped  vector  in  Equation  (3.11)  is 
within  an  error  vector  2zre;*  of  j8q,  where  e;*  is  an  integer  vector  in  K  +  L.  RSC  algorithms  are 

/s  s|e  . 

fundamentally  blind  to  such  errors;  distinct  unwrappings  j 6  and  /30  both  produce  solutions 
to  Equation  (3.5)  in  the  noiseless  case,  as  do  their  respective  projections  onto  K  +  L,  fiK  ( L 
and  /3q  k+l,  in  the  noisy  case.  The  property  of  wrap-invariance  introduced  in  this  chapter 
ensures  that  the  effect  of  such  an  error  on  the  phase  of  the  resulting  RSC  solution  Trsc  is 
either:  merely  a  benign  integer  multiple  of  2n  (c.f.  Section  3.5.2),  or  a  linear  gradient  in  the 
estimated  Fourier  phases,  which  is  equivalent  up  to  an  extra  image  shift  (c.f.  Section  3.7)  in 
the  reconstructed  image.  In  order  to  develop  conditions  for  wrap-invariance  as  they  relate 
to  pattern  design,  we  must  first  characterize  the  effect  of  the  residual  wrap  error  on  the  RSC 

2Without  loss  of  generality,  we  have  selected  e/,  =  0  in  Equation  (3.10)  so  as  to  simplify  the  diagram 
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Noiseless 


Noisy 


P  (wrapped) 


Figure  3.6:  Illustration  of  the  fundamental  ambiguity  of  2n-periodicity  in  RSC  imaging.  Distinct  unwrap¬ 
pings  ft  and  0*  both  produce  solutions  to  Equation  (3.5)  in  the  noiseless  case,  as  do  their  respective  projections 
onto  K  +  L,  PK+L  and  j3*0  K+L,  in  the  noisy  case. 


least-squares  solution.  This  is  the  aim  of  the  next  section. 

3.5.2  Quantifying  the  Effect  of  the  Fundamental  Ambiguity 

Let  us  now  quantify  the  effect  of  this  unresolved  wrap  error  on  the  ultimate  least-squares 
solution,  which  can  be  accomplished  in  two  easy  steps.  Following  standard  least-squares 
principles,  we  first  find  the  projection  fi*K+L  of  the  unwrapped  /3  onto  K  +  L  whose  weighted 
distance  from  j 6  is  minimized.  Note  that  the  e;*  term  in  fi  is  already  in  K  +  L  and  hence 
unchanged  by  projection.  Hence  we  obtain: 


Fk+l  =  W-i(I-Pw)W*j&*=j8( 


'0  ,K+L  +  ^  71  eh 


(3.13) 


where  /Sq  k+l  is  the  projection  of  (I  onto  K  +  L.  We  then  solve  the  system: 


M 


e 

<P 


=  fi. 


K+L 


(3.14) 


As  aforementioned,  M  is  rank-deficient  (by  3),  and  hence  there  will  be  infinitely-many 
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solutions  to  this  system.  Successful  recovery  of  the  Fourier  phases  modulo  2n  requires  a 
solution  preserving  the  integrality  of  the  error  term  ej*.  In  this  case,  we  obtain  a  final  RSC 
solution  27rc  away  from  the  true  solution  for  some  integer  vector  c.  To  achieve  this,  we  rely 
on  the  integer  matrix  decomposition  known  as  the  Smith  Normal  Form,  which  is  described  in 
the  following  Theorem  and  Lemma: 

Theorem  3.5.2.  (Smith  Normal  Form)  (Smith,  1861):  Let  A  be  a  nonzero  m  x  n  integer  matrix 
zvith  rank  r.  There  exist  integer  and  unimodular  3  (and  thus  invertible)  matrices  m  x  m  and  n  x  n 
matrices  U  and  V  respectively  such  that  the  matrix  product  D  =  UAV  is  a  diagonal  matrix  zvhose 
diagonal  entries  D„  (the  so-called  elementary  divisors)  are  non-zero  integers  for  i  <  r,  and  zero  for 
i  >  r.  Moreover,  the  matrices  U  land  respectively,  the  matrix  Vj  represent  the  rozv  land  column I 
operations  which  zeroize  A  below  land  above I  the  diagonal. 

Theorem  3.5.3.  (Elementary  Divisors)  (Smith,  1861):  The  product  of  the  elementary  divisors  is  the 
greatest  common  divisor  (gcd)  of  all  r  x  r  minors  of  A. 

The  proof  of  the  Theorems  above  can  be  found  in  textbooks  such  as  Newman  (1972).  □ 

Let  us  compute  the  Smith  Normal  Form  (SNF)  {U,  D,  V}  of  our  matrix  M.  Let  Um  =  U  1 
and  Vm  =  V  1  so  that  we  can  write: 


M  =  UmDmVm 

where  the  r  diagonal  elements  { d,}  of  Djy  are  the  elementary  divisors  of  M. 
We  can  now  re-write  Equation  (3.14)  above  as: 


(3.15) 


—  fiK+L 

Let  us  choose  the  following  particular  solution  to  Equation  (3.16): 

TrSC  =  Pk+l 

3  A  unimodular  matrix  is  one  whose  determinant  is  ±1 


DmV 


9 

<P 


(3.16) 


(3.17) 
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where  denotes  the  pseudo-inverse  of  D.  The  resulting  error  is  then: 

eRSC  =  V^D+U  J}e*h  (3.18) 

Lemma  3.5.4.  Let  u  =  UM'  e*h  .  The  residual  wrap  error  e^sc  equals  0  mod  2rc  if  and  only  if 
mod  df)  =  0 ,Vz  <  r.  The  proof  of  this  Lemma  is  an  adaptation  of  a  standard  proof  which  can  he 
found  in  most  textbooks  covering  linear  Diophantine  equations  (see,  e.g.  Newman  (1972)). 

□ 

From  this  Lemma,  the  following  Corollary  is  clear: 


Corollary  3.5.5.  (Sufficient  condition  on  SNF  of  RSC  matrix  for  wrap-invariance):  If  the 

elementary  divisors  of  the  measurement  matrix  M  corresponding  to  a  certain  aperture  pattern  are  all 
1,  the  RSC  solution  defined  by  t rsc  is  immune  to  phase-wrapping  error. 

□ 

RSC  patterns  consisting  of  apertures  placed  randomly  on  a  Cartesian  grid  appear  to 
satisfy  this  sufficient  condition  with  high  probability.  We  conducted  a  simple  experiment 
in  which  15  apertures  were  randomly  placed  on  a  10  x  10  grid.  Out  of  256  placements,  66 
were  valid  RSC  patterns  (i.e.  possessed  at  least  critical  redundancy),  and  of  these,  only  2 
had  non-unity  elementary  divisors. 

3.5.3  Relation  to  closure-phase  approaches 

The  SNF  has  been  been  applied  to  the  RSC  phase  problem  before  (Lannes,  2003).  Whereas 
we  have  chosen  to  apply  SNF  directly  to  the  baseline  measurement  matrix,  the  approach 
taken  by  Lannes  (2003)  is  to  instead  treat  the  piston-invariant  phases  of  the  bispectra  (the 
so-called  closure  phases)  as  the  fundamental  observables  from  which  the  object  phases  can  be 
inferred  via  the  relation: 


Co—fC  0  =  yci+2Tzed 


(3.19) 
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where  yc/  are  the  wrapped  closure  phases,  ec/  is  the  corresponding  wrap  vector,  and  C0_j.c  is 
the  matrix  mapping  the  distinct  object  phases  in  the  array  to  closure  phases. 

Using  closure  phases  as  observables  can  be  advantageous  in  low-light  scenarios  in  which 
there  is  not  sufficient  SNR  in  a  single  atmospheric  coherence  time  to  reliably  measure  the 
baseline  phases.  To  overcome  this  low  per-frame  SNR,  atmosphere-invariant  observables 
such  as  the  bispectra  can  be  integrated  over  many  frames  to  build  sufficient  SNR,  and  their 
respective  closure  phases  used  as  reliable  phase  measurements.  Since  the  baseline  phases 
are  known  modulo  2n,  the  linear  combinations  of  them  which  comprise  the  closure  phases 
are  also  known  modulo  2zr. 

Lannes  (2003)  applies  the  SNF  to  the  closure  matrix  C0_>.c.  By  direct  analogy  to  Corollary 
3.5.5,  note  that  if  the  elementary  divisors  of  Coc  are  all  1,  then  the  RSC  solution  is  immune 
to  phase-wrapping  error.  Otherwise,  severe  distortion  is  possible  in  the  resulting  image 
reconstruction.  In  order  to  relate  this  condition  to  Corollary  3.5.5,  let  us  first  define  another 
closure  matrix  Cm^c  which  instead  maps  the  phase  measurements  to  closure  phases.  This 
mapping  consists  of  equations  of  the  form: 

1/123  =  fill  +  @23  ~  fil3  (3.20) 

where  1/123  is  the  closure  phase  associated  with  apertures  1,  2,  and  3,  and  the  f5tj  are  the 
associated  baseline  phases  (see  Equation  (3.2)).  Of  the  (^)  possible  closure  phases,  at  most 
(N21)  can  be  linearly-independent  (see  e.g.  Readhead  et  al.  (1988)).  One  commonly-chosen 
set  of  such  linearly-independent  relations  consists  of  all  the  closure  triangles  involving 
a  given  aperture  A,  and  this  is  the  set  selected  by  Lannes  (2003).  CmH >c  is  therefore  an 
(N2  l)  ><  (2)  matrix.  Lannes  (2003)  accordingly  provides  a  convenient  grouping  of  the 
baselines  into  two  categories:  (1)  spanning  tree  baselines  which  connect  aperture  A  to  all 
other  apertures,  and  (2)  loop  entry  baselines  which  provide  the  closure  for  these  spanning 
tree  baselines.  This  categorization  is  depicted  in  Figure  3.7. 
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Figure  3.7:  Distinction  between  spanning  tree  baselines  (thick,  solid)  and  loop  entry  baselines  (thin,  dotted) 


Given  this  categorization,  we  can  decompose  C „,-+c  into  corresponding  blocks  as: 


C 


m— >c 


Cm— >c  ^(n2~1) 


(3.21) 


where  Cm->c  contains  the  spanning  tree  contributions  to  the  matrix  (which  appear  in 
multiple  closures),  and  I^v  i;  is  the  (V)  x  (V)  identity  matrix  representing  the  loop- 
entry  contributions  (each  of  which  appears  in  only  one  closure).  The  following  property 
follows  from  this  block  form  expression: 


Lemma  3.5.6.  The  elementary  divisors  o/Cm_>.c  are  all  1. 

Proof:  Since  we  have  chosen  a  linearly-independent  subset  of  closure  relations,  r  = 
rank(Cm^c)  =  (N2'1).  There  exists  a  r  x  r  minor  (namely,  I^n-i^)  which  is  equal  to  1. 
Therefore  the  gcd  of  all  r  x  r  minors  is  1,  and  therefore  from  Theorem  3.5.3  (Elementary 
Divisors),  all  elementary  divisors  must  be  1.  □ 

Let  us  now  relate  Cm_>.c  to  the  matrix  C0_>.c  used  by  Lannes  (2003).  Recall  from  the 
discussion  of  bispectra  in  Section  3.4  that  the  closure  relations  eliminate  piston  differences 
in  the  measurements  so  that  C,„_>c  annihilates  the  subspace  L,  i.e.  the  space  spanned  by 
the  columns  of  M  corresponding  to  (p.  Defining  Mg  as  the  submatrix  of  M  containing  the 
columns  corresponding  to  9,  we  have 


0 

0 

Cm— >cM 

= 

<P 

<P 

(3.22) 
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where  C(,_>c  =  Cm_>.cM g.  C 0_>c  is  an  (  N2  1 )  x  d  matrix  which  is  rank-deficient  by  two  4. 

Then  by  direct  analogy  to  Equation  (3.14),  we  can  obtain  valid  RSC  object  phase  solutions 
by  solving: 


0 

r 

0 

Cm^cM 

= 

C»->c  o 

<P 

- 

<p 

=  y*d  +  2ne 


* 

h,cl 


(3.23) 


where  y*cl  is  the  true  unwrapped  closure  vector  and  27rej*c/  is  the  residual  integer  wrapping 
error  vector  after  applying  the  Babai  NP  algorithm  to  solve  the  CVP  problem  associated  with 
matrix  Cm_>.cM0  and  2nec[.  Note  that  if  we  find  a  vector  9*  satisfying  Equation  (3.23),  it  will 
clearly  also  satisfy  the  relation  C0  ^c0*  =  y*;  +  2ne*h  d  of  Lannes  (2003).  Note  furthermore 
that  we  can  solve  the  equation  above  in  two  separate  integer-preserving  steps  of  the  form 
of  Equation  (3.17),  the  first  involving  the  SNF  decomposition  of  Cm_j ,c,  and  the  second 
involving  that  of  M.  Since  the  elementary  divisors  of  Cm_>.c  are  all  1  by  construction  (by 
Lemma  3.5.6)  and  hence  the  first  step  is  thus  integer-preserving,  wrap-invariance  again 
amounts  to  whether  or  not  all  elementary  divisors  of  M  are  1.  Therefore  we  have  the 
following  Proposition  relating  wrap  invariance  for  closure  measurements  to  that  for  raw 
phase  measurements: 


Proposition  3.5.7.  (Sufficient  condition  for  wrap-invariance  of  closure-based  RSC):  If  the 

elementary  divisors  of  the  phase  measurement  matrix  M  are  all  1,  then  the  closure-based  RSC  solution 
zvill  be  wrap-invariant. 


Proof:  (see  Appendix  B.l)  □ 

We  remark  in  passing  that  although  the  preceding  analysis  was  presented  in  the  context 
of  the  traditional  three-aperture  closure,  it  applies  directly  to  the  case  of  closures  involving 
an  arbitrary  number  of  sides.  As  an  example,  consider  the  pattern  shown  in  Figure  3.8.  A 
spanning  tree  for  the  pattern  consisting  of  the  short  baselines  in  an  array  is  depicted.  Let 
{(psp}  denote  the  aperture  phase  differences  in  these  N  —  1  spanning  tree  baselines.  Note 

4To  see  this,  note  that  each  solution  set  to  Equation  (4.49)  above  remains  valid  after  replacing  each  0jj  with 

0-j  =  Oij  ~  z  '  (*>•  -  r i) 
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Figure  3.8:  Bootstrapping  phase  of  a  low-SNR  baseline  (green)  with  subset  (blue)  of  high-SNR  baselines  from 
spanning  tree  baselines  (black) 

that  all  aperture  phase  differences  in  the  array  can  be  expressed  as  linear  combinations 
of  the  {< pSp }.  If  the  aperture  phase  differences  are  known  reliably  via  measurements  of 
the  spanning  tree  baselines,  we  can  use  these  measurements  to  cancel  the  aperture  phase 
differences  in  all  other  measurements  (of  which  one  example  is  shown  in  green).  The 
idea  of  using  such  generalized  closure  phases  (Martinache,  2010)  is  indeed  the  mathematical 
foundation  for  the  promising  technique  known  as  baseline  bootstrapping  in  which  high-fidelity 
phase  measurements  of  several  high-SNR  baselines  (typically  the  short  baselines)  to  cancel 
the  atmosphere  on  each  lower-SNR  (and  hence  lower  fidelity)  baseline. 

Note  that  for  an  arbitrary  N-aperture  pattern,  there  will  in  general  be  N  —  1  spanning 
tree  baselines  and  thus  (^)  —  (N  —  1)  =  ( ,V2  1 )  generalized  closures,  each  involving  a 
distinct  closing  (or  loop-entry)  baseline.  Therefore  the  resulting  measurement  matrix  can  be 
expressed  exactly  as  in  Equation  (3.21)  above  and  hence  the  preceding  analysis  holds. 

While  this  section  has  considered  a  few  possibilities  for  phase  observables,  relating 
mathematical  conditions  for  wrap  invariance  to  a  physical  condition  on  aperture  placement 
is  more  intuitive  when  considering  the  raw  phase  measurements  as  opposed  to  their  closures. 
For  this  reason,  for  the  remainder  of  the  chapter  we  will  work  directly  with  the  baseline 
phases  and  their  associated  wrapping  errors.  In  particular,  we  will  begin  by  connecting 
these  wrapping  errors  with  analogous  ambiguities  in  recently-developed  phasor-based 
approaches. 
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3.6  Implications  of  Wrap  Ambiguities  on  Pattern  Design 


In  this  section  we  use  the  mathematically-sufficient  conditions  for  wrap  invariance  from 
the  previous  section  to  show  that  aperture  patterns  whose  interferometric  graph  satisfies  a 
certain  loop-free  condition  are  wrap-invariant.  Here  we  define  the  interferometric  graph  in 
the  standard  way:  it  is  simply  the  graph  formed  by  connecting  the  array's  apertures  (the 
nodes  of  the  graph)  with  edges  representing  the  array's  baselines. 

This  condition  is  founded  on  the  Theorem  3.5.3  (Elementary  Divisors)  in  Section  3.5  and 
the  following  definition  of  the  matrix  determinant.  This  definition  is  given  in  many  linear 
algebra  texts  (see  e.g.  Bretscher  (2001)). 

Definition  3.6.1.  Matrix  Determinant  in  terms  of  patterns:  Suppose  we  have  an  n  x  n 
matrix  A.  Define  a  pattern  as  a  selection  of  n  entries  of  the  matrix  in  which  there  is  only  one 
chosen  entry  in  each  row  and  one  in  each  column  of  the  matrix.  Furthermore,  we  denote  a 
pair  of  numbers  in  a  pattern  as  inverted  if  one  of  them  is  located  above  and  to  the  right  of 
the  other.  Then  we  can  obtain  the  determinant  of  A  by  summing  the  products  associated 
with  all  patterns  with  an  even  number  of  inversions  and  subtracting  the  products  associated 
with  all  the  patterns  with  an  odd  number  of  inversions. 

Consider  a  r  x  r  sub-matrix  Mj  consisting  of  a  set  I  of  linearly-independent  rows  in  M 
and  d  +  N  —  3  linearly-independent  columns.  Our  goal  will  be  to  find  conditions  under 
which  such  a  sub-matrix  contains  only  one  pattern  with  a  non-zero  product,  in  which  case 
the  determinant  will  be  ±1  from  the  definition  above.  This  will  guarantee,  by  the  Theorem 
3.5.3  (Elementary  Divisors)  of  the  previous  section,  that  the  elementary  divisors  of  M  are  all 
1,  which  will  in  turn  ensure  that  the  error  in  the  RSC  solution  Trsc  will  be  0  mod  2tc  by 
Corollary  3.5.5. 

We  begin  our  search  for  such  a  matrix  by  forming  the  intermediate  r  x  (r  +  3)  matrix 
M/  consisting  of  I  linearly-independent  rows  in  M,  and  define  the  corresponding  sub-graph 
G  containing  only  the  baselines  corresponding  to  these  rows.  In  particular,  let  us  examine 
the  pattern  restrictions  encountered  if  we  select  a  given  column  in  M/  for  participation  in 
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our  r  x  r  sub-matrix  Mj.  Namely,  we  will  sequentially  identify  those  special  matrix  entries 
which  must  be  part  of  a  pattern  with  a  non-zero  product.  For  example,  all  non-redundant 
measurements  contain  a  singleton  ±1  in  the  column  associated  with  their  object  phase.  All 
non-zero  patterns  must  clearly  contain  this  ±1  and  so  we  can  select  these  singleton  object- 
phase  entries  as  guaranteed  participants  in  a  non-zero  pattern.  Moreover,  all  measurements 
containing  a  leaf  node  (i.e.  a  node  with  a  single  connection)  in  G  contain  a  singleton  ±1  in 
the  column  associated  with  their  leaf  node.  All  non-zero  patterns  must  clearly  contain  this 
±1  as  well.  Thus  we  can  also  select  these  leaf  node  entries  as  guaranteed  participants  in  a 
non-zero  pattern. 

There  may  be  cascading  implications  of  such  singleton  measurements.  To  illustrate  this, 
consider  the  scenario  shown  in  Figure  3.9.  A  simple  RSC  array  is  shown  on  the  left.  A 
subset  of  the  baselines  in  one  possible  linearly-independent  sub-matrix  M;  is  depicted.  A 
simplified  depiction  of  the  matrix  Mj  is  shown  in  which  all  non-zero  entries  have  been 
colored  black  and  all  zero  entries  have  been  colored  white  for  simplicity.  In  Step  A  of  the 
reduction  process,  object  phase  %  is  selected  for  participation  since  it  is  a  singleton  object 
phase.  Its  corresponding  row  (i.e.  row  5)  in  M/  is  then  eliminated  from  participation  since 
the  remaining  entries  in  this  row  cannot  participate  in  a  pattern  (by  definition  of  a  pattern). 
In  Step  B,  the  aperture  6  entry  (p(,  in  row  9  is  selected  since  it  has  become  a  leaf  node  in 
the  pattern,  and  row  9  can  then  be  eliminated.  Then  in  Step  C,  object  phase  04  is  then 
selected  by  virtue  of  becoming  a  singleton  object  phase,  and  row  4  is  then  eliminated.  This 
selection/ elimination  process  can  be  repeated  beyond  the  steps  shown  in  the  Figure,  until 
either  no  leaf  nodes  and  singleton  object  phases  remain,  or  there  are  no  baselines  left  to 
eliminate.  We  formalize  the  pattern  reduction  process  in  the  listing  Algorithm  2. 

We  can  see  that  any  baseline  in  an  interferometric  graph  which  does  not  belong  to  a  loop 
will  be  eliminated  in  the  reduction  process.  The  only  structures  in  the  graph  that  persist 
after  this  reduction  are  sets  of  loops  with  a  certain  property.  Namely  we  define  a  persistent 
loop  set  as  a  set  of  loops  that  contains  at  least  two  instances  of  every  baseline  contained  in 
the  set.  (A  set  can  consist  of  any  number  of  loops,  including  one).  With  this  definition. 
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Figure  3.9:  Example:  Reducing  an  aperture  pattern  and  associated  matrix  to  identify  Persistent  Loop(s) 


Algorithm  2  Pattern  Reduction  Algorithm 

Require:  R  {where  R  is  the  set  of  the  baselines  corresponding  to  M/,  where  I  denotes  the 
indices  of  a  linearly-independent  subset  of  d  +  N  —  3  rows  of  M) 

while  |  R  |  >  0  do 

1.  Leaf  Nodes 

1.1  remove  any  remaining  baselines  containing  leaf  nodes  from  R 

1.2  add  the  associated  apertures  to  the  list  N 

2.  Singleton  Object  Phases 

2.1  remove  any  remaining  baselines  containing  singleton  object  phases  from  R 

2.2  add  the  associated  object  phases  to  the  list  O 

if  no  baselines  removed  in  the  current  iteration  then 
return  PERSISTENT 

end  if 
end  while 
return  LOOPFREE 
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absolute  invariance  may  be  possible  if  the  graph  of  the  redundant  baselines  does  not  contain 
any  persistent  loop  sets.  Algorithm  2  returns  PERSISTENT  if  persistent  loops  exist  and 
LOOPFREE  if  the  pattern  is  completely  reduced  and  therefore  free  of  persistent  loops. 

Note  that  in  the  latter  case,  we  will  have  eliminated  r  rows  from  M;.  Since  each  baseline 
elimination  is  associated  with  object  phase  or  aperture  selection  from  distinct  columns,  and 
Mj  contains  r  +  3  columns,  there  will  be  exactly  three  extraneous  columns  not  involved 
in  the  reduction  process.  The  r  x  r  submatrix  obtained  by  selecting  the  non-extraneous 
columns  (i.e.  those  corresponding  to  the  selected  object  phases  and  leaf  nodes  in  O  and 
N,  respectively)  will  then  have  a  determinant  of  ±1  by  virtue  of  having  a  single  non-zero 
pattern  revealed  by  the  reduction  process.  Having  ensured  the  existence  of  a  unit  r  x  r 
minor,  the  gcd  of  all  r  x  r  minors  must  be  1.  We  have  hence  confirmed  the  elementary 
divisors  must  be  all  l,  and  thereby  ensured  that  the  pattern  is  wrap-invariant  (c.f.  Theorem 
3.5.3  (Elementary  Divisors),  and  Corollary  3.5.5).  We  summarize  the  sufficient  condition  as 
follows: 

Proposition  3.6.1.  (Sufficient  conditions  on  aperture  pattern  for  wrap-invariance):  Consider 
the  graph  of  an  aperture  pattern  which  contains  d  distinct  baselines  and  any  set  ofN  —  3  linearly- 
independent  redundant  baselines.  If  this  graph  does  not  contain  persistent  loop  sets  (in  the  sense 
defined  above),  the  matrix  Mj  formed  by  these  independent  measurements  zvill  have  determinant  ±1. 
As  a  result  Corollary  3.5.5  zvill  hold,  thereby  guaranteeing  that  RSC  solution  Trsc  be  invariant 
to  wrapping  of  the  phase  measurements. 

We  have  hence  arrived  at  a  a  physical  definition  of  a  zvrap -invariant  pattern.  We  now 
apply  Algorithm  2  to  the  example  pattern  shown  in  Figure  3.2.  Algorithm  2  reduces  the 
pattern  to  the  persistent  loop  set  shown  in  Figure  3.10. 

The  elementary  divisors  of  the  pattern's  measurement  matrix  are  not  all  1;  they  are  all  1 
except  for  a  singleton  3  and  hence  det( Mj)  mod  3  =  0  for  all  choices  of  M/. 

Having  traced  the  distortion  induced  by  phase  wrapping  to  physical  property  of  the 
array  itself,  we  now  return  to  Figure  3.3.  The  lower  left  panel  shows  the  reconstruction  result 
with  the  SNF-based  Phase  method  described  in  Section  3.5.  The  closure  phase  approach 
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Persistent  Loop 


Figure  3.10:  Persistent  Loop  set  at  center  of  pattern  in  Figure  3.2 

yields  the  same  corruption  in  reconstruction,  as  the  elementary  divisors  of  Coc  are  also  all  1 
except  for  a  singleton  3. 

There  are  several  simple  ways  to  amend  this  pattern  so  that  it  is  wrap-invariant.  While 
the  most  intuitive  of  these  involve  moving  the  apertures  involved  in  the  persistent  loop 
shown  in  Figure  3.10,  these  approaches  leave  gaps  in  the  UV-sampling  pattern.  An  alternate 
approach  that  preserves  the  UV-sampling  is  to  add  an  aperture  to  the  center  of  the  pattern 
as  shown  in  Figure  3.11.  This  results  in  additional  linear ly-independent  redundancies 
colored  in  blue  and  green,  respectively,  in  the  Figure.  These  additions  replace  baselines  in 
the  persistent  loop,  allowing  this  loop  to  be  broken.  With  wrap-invariance,  reconstruction 
results  match  the  true  image  in  both  the  phase  and  phasor  approaches  as  respectively  shown 
in  Figure  3.12.  In  the  top  row,  reconstruction  results  are  displayed  for  the  phase  (left)  and 
phasor  (right)  approaches  for  the  noiseless  case.  The  image  distortion  present  in  Figure 
3.3  has  been  completely  eliminated  by  tweaking  the  pattern  so  that  it  is  wrap-invariant. 
Analogous  results  for  an  SNR  of  25  dB  are  displayed  in  the  bottom  row.  Flere  we  define 
SNR  as  the  ratio  of  the  phasor  magnitude  at  visibility  1  (i.e.  zero  spatial  frequency)  to  the 
standard  deviation  of  the  noise,  which  we  have  assumed  to  be  complex  Gaussian  and  i.i.d. 
across  spatial  frequency  for  this  simulation. 
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Figure  3.11:  Amended  Pattern 


Figure  3.12:  Reconstruction  Results  for  Amended  Pattern 
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3.7  Wrap-invariance  and  Practical  RSC  Calibration 


In  previous  sections  we  established  that  the  Phase  Approach  can  be  made  robust  to  phase¬ 
wrapping  using  the  Smith  Normal  Form  (SNF)  and  algorithms  from  lattice  theory.  Moreover, 
the  SNF  provided  a  mathematical  framework  for  the  notion  of  wrap-invariance.  From  a 
practical  standpoint,  however,  computation  of  the  Smith  Normal  Form  is  likely  to  become 
a  computational  burden  for  large  arrays,  as  in  those  under  current  consideration  with 
Nap  ~  102  and  n  ~  104  baselines  (Zheng  et  al.,  2014).  Techniques  not  requiring  such  a 
computation  are  hence  of  strong  practical  interest.  In  this  section  we  show  that  the  wrap- 
invariance  property  checked  by  Algorithm  2  provides  a  certificate  for  reliable  reconstruction 
with  such  techniques,  in  the  presence  of  the  fundamental  271-periodicity  of  interferometric 
measurements. 

3.7.1  Practical  Phase  Approaches 

Our  approach  here  will  rely  upon  the  well-known  Singular  Value  Decomposition  (SVD)  of 
the  measurement  matrix  M,  which  is  given  by: 

M  =  (3.24) 

in  which  U,r  and  V,r  are  m  x  rn  and  (d  -\-  N)  x  (d  +  N  )  orthogonal  matrices,  respectively. 
L,r  is  a  in  x  (d  +  N)  diagonal  matrix  with  r  non-zero  diagonal  entries  (the  so-called  singular 
values  of  M),  where  r  =  rank(M)  =  d  +  N  —  3. 

Lemma  3.7.1.  The  final  3  columns  of  Va  form  a  basis  for  the  nullspace  of  M. 

Proof:  This  follows  from  the  fact  M  is  rank-deficient  by  3,  and  standard  properties  of  the 
right  singular  vectors  comprising  in  the  SVD.  (Bretscher,  2001)  □ 

Now  recall  that  in  Section  3.4,  we  provided  an  one  particular,  SNF-based,  solution  to 
Equation  (3.14).  Flere  we  instead  consider  the  complete  set  of  solutions  to  this  Equation, 
which  is  given  by: 
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(3.25) 


TV  =  M  +  CP*o,k+l  +  2/r  e*h)  +  r0 

where  Tq  is  any  vector  in  the  nullspace  of  M,  and  M+  denotes  the  pseudo-inverse  of  M, 
whose  matrix  elements  are  derived  directly  from  the  SVD  above. 

M+  =  V^E+Uj  (3.26) 

where  E+  is  a  (d  +  N)  x  m  diagonal  matrix  whose  r  non-zero  diagonal  entries  are  the 
reciprocals  of  the  corresponding  non-zero  entries  in  E(r. 

Typically  Phase  Approach  techniques  implicitly  select  a  particular  solution  from  the 
family  of  solutions  in  Equation  (3.25)  via  augmenting  the  matrix  M  with  additional  con¬ 
straints.  The  most  common  among  these  (Wieringa,  1992)  (Wijnholds  and  Noorishad,  2012) 
is  to  enforce  that:  E  <P  =  0,  <pt r,-  =  0,  where  r,  is  the  vector  position  of  the  z-th  aperture  in 
the  array. 

Note  that  the  error  resulting  from  application  of  this  pseudo-inverse  to  the  unwrapped 
measurement  vector  will  be  given  by: 

Ine^  =  M+(27ie;*)  (3.27) 

Note  furthermore  that  the  vector  representing  the  solution  error  ea  has  two  parts:  one 
for  the  error  in  the  Fourier  phases  (which  we  denote  with  the  index  K),  and  one  for  the 
error  in  the  atmosphere /piston  (which  we  will  denote  with  the  index  L),  i.e.: 

eU  =  [ecr,Kr  (3.28) 

We  will  focus  attention  on  the  latter,  since  it  is  of  direct  relevance  for  image  formation. 
Let  us  express  the  spatial  frequencies  measured  by  an  array  as  two-element  vectors  of  the 
form  (ojx,  LOy  ).  Let  X  be  the  d  x  2  matrix  containing  these  spatial  frequencies.  Note  then  that 
the  phase-wrap  error  will  manifest  itself  merely  as  an  image  shift  if  and  only  if  this  error  is 
a  (modulo-27r)  phase  ramp,  i.e.  there  exists  a  2-element  shift  vector  z  and  an  integer  vector 
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k  which  together  satisfy: 


2nea,K  ~  2/rXz  =  Ink  (3.29) 

Substituting  from  Equation  (3.27)  we  obtain: 

M+  (2/re)  -  2/rXz  =  Ink  (3.30) 

where  M(<:  denotes  the  sub-matrix  of  M+  formed  by  the  rows  associated  with  K. 

Dividing  through  by  2n  we  obtain  the  equation:  M(<  e/(  —  Xz  =  k.  Note  that  each 
element  of  can  be  expressed  as  some  rational  number  h.  Similarly  we  first  assume 
X  contains  rational  spatial  frequencies  with  greatest  common  denominator  qx.  Then  we 
can  multiply  through  by  the  least-common-multiple  (LCM)  of  the  {qi}  and  qx  to  obtain  a 
system  of  equations  whose  coefficients  are  guaranteed  to  be  integer  (i.e.,  we  have  a  linear 
Diophantine  system).  Let  this  LCM  be  denoted  as  l.  Then  we  have,  after  rearranging  terms, 

/Xz  =  l(M+e*h  -k)  (3.31) 

We  now  wish  to  determine  conditions  under  which  there  exist  vectors  k  and  z  satisfying 
this  overdetermined  Diophantine  system.  Applying  the  Smith  Normal  Lorm  decomposition 
(c.f.  Theorem  3.5.2)  to  the  matrix  IX  this  time,  and  noting  that  rank(X)  =  2,  we  have: 


Dx  =  Ux(/X)Vx  (3.32) 

where  Ux  and  Vx  are  unimodular  matrices  of  size  d  x  d  and  2x2,  respectively,  and  Dx 
is  a  rectangular  diagonal  matrix  whose  entries  are  zero  below  row  2. 

If  we  left-multiply  Equation  (3.31)  by  Ux  on  both  sides,  we  obtain: 

!UxXz  =  /Ux(M+e;:  -  k)  (3.33) 

Using  Equation  (3.32)  and  the  fact  that  Vx  is  a  unimodular  (and  hence  invertible)  matrix, 
we  can  then  write: 
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DxVx1z  =  /(UxM+e;i-Uxk)  (3.34) 

We  are  now  in  position  to  prove  the  main  result  of  this  section,  which  is  preceded  by  the 
following  Lemma: 

Lemma  3.7.2.  Given  wrap-invariance,  the  vector  UxMxe;*  has  integer  entries  below  row  2. 

Proof:  (see  Appendix  B)  □ 

Proposition  3.7.3.  If  a  pattern  is  wrap-invariant  (in  the  sense  of  Section  3.5),  reconstruction  error 
induced  by  phase  wrapping  is  limited  to  an  image  shift. 

Proof:  We  re-arrange  Equation  (3.34)  above  so  that  it  reads: 

yDxVx'z  -  GxM+e*h  =  — Uxk  (3.35) 

Let  v  =  yDxV^z  —  UxM^eJ.  Note  that  since  Dx  is  zero  below  row  2,  the  entries  of  v 

below  row  2  will  be  equal  to  those  of  (-UxM^eJj),  which  are  integers  by  Lemma  3.7.2.  Now 
consider  the  first  and  second  entries  of  v.  Let  f  be  the  vector  containing  the  fractional  parts 
of  the  first  two  elements  of  vector  UxM^e(),  and  let  A  be  the  invertible  matrix  consisting 
of  the  first  two  rows  of  }DxVx  1 .  Without  loss  of  generality,  choose  z*  =  A  1  f  so  that  the 
fractional  part  f  is  annihilated,  leaving  only  integer  elements  in  the  first  two  entries  of  v. 
Hence  we  now  have: 


v  =  — Uxk  (3.36) 

with  v  ensured  to  contain  only  integer  elements.  Since  Ux  is  unimodular,  the  vector 
k  =  — Ux  1  v  will  be  integral.  We  have  thus  found  a  pair  (z*,k*)  with  integer  k  which 
satisfies  the  Equation  (3.34).  Since  Equation  (3.34)  is  related  to  Equation  (3.31)  via  a 
unimodular  (and  hence  invertible)  mapping  Ux,  invariance  is  hence  proven.  □ 

With  the  previous  result,  we  have  characterized  the  complete  family  of  Phase  Approach 
solutions  given  in  Equation  (3.25).  Namely  we  have  shown  that,  for  a  wrap-invariant  pattern. 
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the  family  differs  by  at  most  an  image  shift  from  the  true  solution.  Returning  to  our  running 
example  in  Figure  3.11,  we  verified  that  different  solution  choices  from  the  family  given  in 
Equation  (3.25)  simply  resulted  in  shifts  of  an  otherwise  pristine  image  in  the  reconstruction. 
On  the  other  hand  with  wrap-variant  pattern  in  Figure  3.2,  image  distortion  of  the  severity 
of  Figure  3.3  was  again  observed,  as  expected. 

3.7.2  Practical  Phasor  Approaches 

Though  traditional  treatments  employ  the  phase  approach  of  the  previous  section  which 
operates  on  baseline  phases,  recent  papers  (e.g.  Marthi  and  Chengalur  (2014),  Fiu  et  al. 
(2010))  have  shown  that  approaches  which  operate  at  the  phasor  level  can  be  superior 
in  accuracy.  Fiu  et  al.  (2010)  developed  a  Gauss-Newton-type  Non-linear  Feast-Squares 
(NFS)  solver  and  showed  it  produced  unbiased  phase  estimates,  in  contrast  with  the 
biased  ones  provided  by  the  phase  approach.  Marthi  and  Chengalur  (2014)  and  Wijnholds 
and  Noorishad  (2012)  have  also  proposed  low-complexity  phasor-based  approaches  and 
demonstrated  performance  near  the  Cramer-Rao  Bound.  Though  the  capacity  of  the  Phasor 
approaches  to  produce  superior  accuracy  relative  to  Phase  approach  has  been  demonstrated, 
the  former's  convergence  issues  can  be  mitigated  via  initialization  with  the  results  of  the 
latter  (see  e.g.  Fiu  et  al.  (2010),  and  Zheng  et  al.  (2014)). 

The  implementations  of  the  Phasor  Approach  typically  employ  the  following  measure¬ 
ment  model: 


Vij  =  gig]fij  +  nij  (3.37) 

where  Vtj  is  the  complex  visibility  observed  between  apertures  i  and  /,  gj  =  \g1\el'l>‘  and 
gj  =  \gj\e^i  are  the  complex  gains  of  these  apertures,  /,,  is  the  true  complex  visibility 
measured  by  this  pair,  and  n  is  complex  measurement  noise.  Note  that  the  phase  difference 
between  g,  and  gj  is  simply  the  optical  path  difference  between  apertures  i  and  j  introduced 
in  the  previous  section.  Given  this  model,  NFS  approaches  attempt  to  find  a  set  of  complex 
phasors  {g,}  and  {/}  which  minimize  an  objective  function  of  the  form: 
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(3.38) 


a  =  LLwMVii-Sigjfii)(yij-shifij)\\ 

i  j>i 

Minimization  of  A  with  respect  to  the  unknowns  (i.e.  distinct  object  and  antenna 
complex  gains)  can  be  accomplished  with  iterative  application  of  the  following  updates, 
as  reported  by  Marthi  and  Chengalur  (2014)  in  the  context  of  radio  interferometry,  and  by 
Lacour  et  al.  (2007)  in  the  context  of  optical  interferometry: 


8k  = 


fb  ~ 


zvkjgjfkjVkj 

YLj^kwkj\gj\  \fkj\ 
Lj>k8*k8jvkj 


(3.39) 


(3.40) 


Lj>kWkj\gk\  \gj\ 

where  the  {/),}  are  the  true  complex  visibilities  of  the  distinct  object  phases  in  the  array. 

Due  to  the  circularity  of  these  definitions,  these  equations  must  be  solved  iteratively. 
Starting  from  an  initial  guess  for  all  phasors.  Equation  (3.39)  is  solved  to  obtain  a  better 
estimate  for  the  {y^}  and  then  these  {yy}  are  used  to  obtain  refined  estimates  of  the  { /), } 
through  Equation  (3.40).  In  the  next  iteration,  these  {//,}  are  used  to  further  refine  {gk},  and 
so  on. 

Though  there  are  other  means  of  minimizing  objectives  of  the  form  A  (Wijnholds  and 
Noorishad,  2012)  (Liu  et  al.,  2010),  we  omit  discussion  of  them  here;  our  present  purpose 
is  to  characterize  the  correctness  of  the  solutions  themselves,  regardless  of  how  they  are 
obtained.  As  has  been  noted  before  (see  e.g.  Lannes  and  Anterrieu  (1999)),  there  are  strong 
connections  between  phase-  and  phasor-based  approaches.  To  see  this,  let  z  be  the  vector  of 
products  {gig*f\i-j\}  which  minimize  A.  We  rewrite  Equation  (3.38)  as: 


A  =  E  E  ™ii  1 1  %  -  -  zij)  1 1  (3-41) 

i  j>i 

Define  r,/  =  e,lntli>Zij  for  an  arbitrary  integer  and  r  as  the  vector  containing  the  r,;. 
Note  that  r  also  minimizes  A  since  the  rotations  {el2nn'i}  do  not  change  the  values  of 
the  residuals  in  A.  Hence  any  set  of  rotated  phasors  {y,}  and  {f^  ;  }  whose  products 
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produce  the  vector  r  will  also  minimize  A.  Note  that  the  set  of  such  valid  phase  vectors 
(i.e.  the  concatenations  of  possible  { Zg, }  and  { Z/|,-_  /| })  includes  the  complete  family  of 
Phase-approach  solutions  t(T  in  Section  3.7.1  with  f"K+  L  =  Zr  (where  Zr  is  the  vector  of  the 
phases  of  the  complex  vector  r).  In  other  words,  the  valid  phase  component  of  the  phasor 
approach  solutions  is  not  unique,  and  the  minimization  of  A  admits  the  same  solution 
ambiguity  depicted  in  Figure  3.6.  Hence  we  see  that  integer  ambiguities  present  in  the 
phase  approach  do  not  disappear  in  the  phasor  approach;  in  fact,  the  unwrapped  candidate 
solutions  of  the  phase-based  approach  correspond  to  minima  of  the  phasor-based  objective. 
In  practice  Phasor  Approach  techniques  may  converge  to  any  one  of  these  minima.  Hence 
the  critical  issue  for  reliable  image  reconstruction  is  again  the  nature  of  the  difference 
between  these  valid  minima  and  the  true  solutions.  Based  on  the  connections  we  have 
drawn  with  Phase-approach  solutions  above,  the  following  Proposition  is  clear: 

Proposition  3.7.4.  If  a  pattern  is  wrap-invariant,  the  global  minima  of  the  Phasor-Approach  objective 
caused  by  the  inherent  Irc-periodicity  in  the  objective's  residual  differ  from  the  true  solution  merely 
by  an  image  shift. 

In  practice  we  see  that,  as  in  the  Phase  approach,  if  the  pattern  is  not  wrap-invariant, 
the  Phasor  Approach  suffers  from  false  global  minima  producing  severe  distortion  in 
the  resulting  reconstruction.  The  lower  right  panel  of  Figure  3.3  shows  the  resulting 
reconstruction  produced  by  the  Phasor  method  using  the  updates  in  Equations  (3.39)  and 
(3.40).  To  show  the  correspondence  in  solutions  between  two  types  of  approaches,  we 
provided  the  result  of  the  Phase  approach  as  the  initial  point  for  the  Phasor  approach 
as  is  common  in  practice  (Liu  et  al.,  2010),  (Zheng  et  al.,  2014).  Indeed  this  point  is  a 
global  minimum  of  A  which  the  updates  above  cannot  escape,  and  as  a  result  we  observe 
virtually-identical  distortion  to  that  of  the  Phase  approach. 

We  repeated  the  experiment  with  the  wrap-invariant  pattern  in  Figure  3.11,  and  as 
expected,  the  results  were  pristine  in  the  noiseless  case  and  virtually-identical  to  those  of 
the  Phase  approach  in  the  noisy  case.  For  completeness  the  results  are  given  in  Figure  3.13. 

As  expected,  different  initializations  of  the  updates  in  Equations  (3.39)  and  (3.40)  simply 
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Phasor  Approach  (Noiseless)  Phasor  Approach  (SNR  =  25  dB) 


Figure  3.13:  Reconstruction  Results  for  Phasor  Approach 


resulted  in  shifts  of  an  otherwise  pristine  image  in  the  reconstruction. 

3.8  Conclusions 

In  this  chapter,  we  have  examined  the  ambiguities  caused  by  the  2/r-periodicity  of  inter¬ 
ferometric  phase  in  Redundant  Spacing  Calibration.  In  particular  we  have  described  their 
fundamental  presence  in  existing  RSC  methods  whether  the  observables  considered  are  the 
measured  baseline  phasors  or  their  phases.  In  the  former,  e.g.  by  Greenaway  (1990),  they 
are  manifested  as  phase- wrapping  errors,  and  in  the  latter,  e.g.  by  Marthi  and  Chengalur 
(2014),  as  false  minima.  We  have  demonstrated  that  in  either  case,  these  ambiguities  can 
result  in  noticeable  distortion  of  the  reconstructed  image. 

Using  the  Closest- Vector-Problem  formulation  of  the  unwrapping  problem  due  to  Lannes 
and  Anterrieu  (1999),  we  have  developed  the  notion  of  a  wrap-invariant  pattern.  For  wrap- 
invariant  patterns,  the  impact  of  the  2zr-periodicity  can  be  completely  eliminated  (c.f.  Section 
3.5)  using  well-known  algorithms  from  lattice  theory  and  the  Smith  Normal  Form,  and 
reduced  to  a  mere  image  shift  (c.f.  Section  3.7)  when  existing,  fast  approaches  are  used. 
Phase-approach  solutions  (Arnot  et  al,  1985),  (Greenaway,  1990),  (Wieringa,  1992),  (Lannes 
and  Anterrieu,  1999))  are  commonly  used  to  quickly  obtain  an  initial  point  to  aid  in  the 
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convergence  of  the  Phasor  approach.  They  are  obtained  by  selecting  specific  solutions 
from  the  general  family  in  Equation  (3.25)  via  enforcement  of  additional  constraints  on 
the  solution.  Phasor  approaches  seek  those  complex  gains  and  object  visibilities  which 
minimize  a  squared-residual  with  respect  to  the  observed  complex  visibilities  (Wijnholds 
and  Noorishad,  2012),  (Liu  et  al.,  2010),  (Marthi  and  Chengalur,  2014)  using,  for  example, 
gradient-descent  methods.  We  have  seen  that  the  same  27r-ambiguity  which  creates  a  family 
of  Phase-approach  solutions  also  produces  a  corresponding  family  of  global  minima  in  the 
Phasor  approach.  For  either  case,  our  results  show  that  for  wrap-invariant  patterns,  this 
family  represents  merely  shifted  versions  of  the  true  image.  We  have  also  extended  this 
analysis  to  show  that  wrap-invariant  patterns  admit  reliable  imaging  using  standard,  and 
generalized,  phase  closures.  Conversely  we  show  that  patterns  which  are  not  wrap-invariant 
can  suffer  from  distortion  of  the  sort  depicted  in  Figure  3.3. 

The  prognosis  for  mitigation  of  the  ambiguity  issues  raised  in  this  chapter  is  quite 
positive.  Random  patterns  appear  to  satisfy  the  wrap-invariance  condition  with  high 
probability.  Moreover,  failure  to  meet  this  condition  amounts  to  the  existence  of  a  particular 
kind  of  cycle  in  the  interferometric  graph  which  can  be  easily  isolated;  a  chief  contribution 
of  this  chapter  is  a  simple  algorithm  for  identifying  such  cycles  so  that  they  can  be  removed 
by  the  array  designer.  Finally,  we  have  shown  an  example  execution  of  the  algorithm  to 
diagnose  a  member  of  a  popular  pattern  which  is  not  wrap-invariant.  It  is  clear  that  with 
careful  array  design,  both  Phase-  and  Phasor-based  RSC  techniques  can  reliably  produce 
quality  image  reconstructions  free  from  discernible  artifacts. 
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Chapter  4 


Robust  image  reconstruction  with 
redundant  arrays  and  generalized 
closure  phases 

4.1  Chapter  Overview 

Redundant  Spacing  Calibration  (RSC)  techniques  employ  redundancy  in  the  baselines  of  a 
telescope  array  to  eliminate  the  contribution  of  atmospheric  turbulence  in  the  interferometric 
observables.  Whereas  conventional  techniques  for  this  elimination  require  the  enforcement 
of  prior  constraints  on  the  underlying  image,  RSC  algorithms  can  be,  in  principle,  mathemat¬ 
ically  well-posed  and  hence  require  no  such  prior  knowledge.  Traditionally  these  algorithms 
have  been  applied  directly  to  the  fringe  measurements.  However,  in  scenarios  of  low  pho¬ 
ton  flux,  such  as  those  arising  in  the  observation  of  dim  objects  in  space,  single-exposure 
fringe  measurements  are  not  reliable  observables  in  general.  Instead  one  must  rely  on 
time-averaged,  atmosphere-invariant  quantities  such  as  the  bispectrum.  In  this  chapter,  we 
develop  a  novel  algorithm  for  redundant  arrays  which  provides  robust  image  reconstruction 
using  integrable  atmosphere-invariant  observables.  Our  algorithm  utilizes  standard  linear 
estimation  methods,  as  well  as  techniques  from  lattice  theory,  to  reliably  estimate  the  Fourier 
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phase.  Moreover,  we  provide  theoretical  and  empirical  evidence  that  generalizing  the  classi¬ 
cal  bispectrum  to  higher-order  atmosphere-invariant  observables,  which  we  call  n-spectra, 
can  offer  significant  performance  gains.  Our  selection  of  an  independent  and  high-SNR  set 
of  n-spectra  leverages  the  notion  of  the  minimum  cycle  basis  from  graph  theory.  We  analyze 
the  expected  shot-noise-limited  performance  of  our  algorithm  for  both  pairwise  and  Fizeau 
interferometric  architectures,  and  corroborate  this  analysis  with  simulation  results  showing 
performance  near  the  Cramer-Rao  bound.  Lastly,  we  apply  techniques  from  the  field  of 
sparse  recovery  (SR)  to  perform  image  reconstruction  from  the  estimated  complex  visibilities. 

4.2  Problem  Statement  and  Prior  Work 

Recall  from  previous  Chapters  that  interferometric  measurements  take  the  form: 


p  =  6  +  A$  (4.1) 

where  A  is  an  (^)  x  N  matrix  whose  rows  compute  the  piston  differences  involved  in  each 
measurement. 

In  Chapter  2  we  introduced  a  class  of  algorithms  which  decouple  the  Fourier  phases  6 
from  the  atmosphere  contributions  (p  by  imposing  prior  constraints.  With  estimates  for  sparse 
samples  of  the  Fourier  Transform,  we  showed  that  the  residual  ill-posed  inverse  problem 
could  be  regularized  effectively  using  techniques  from  sparse  recovery.  As  interferometric 
systems  are  tasked  with  the  observation  of  increasingly-complex  scenes,  approaches  which 
minimize  reliance  on  specialized  prior  assumptions  are  expected  to  gain  appeal.  In  Chapter 
3  we  described  the  family  of  RSC  techniques  as  an  intrinsically  well-posed  alternative  to 
prior-regularized  phase  recovery  is  to  use  baseline  redundancy  to  explicitly  solve  for  piston 
variation.  We  also  examined  the  important  issue  of  integer  phase  ambiguities  in  the  context 
of  these  techniques.  Array  design  considerations  based  on  the  notion  of  wrap-invariance 
(Kurien  et  al.,  2016),  as  well  as  algorithms  from  lattice  theory  (Lannes  and  Anterrieu,  1999), 
can  be  used  to  solve  these  ambiguity  issues. 
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In  this  chapter,  we  leverage  elements  from  both  Chapters  2  and  3  to  develop  a  novel 
algorithmic  framework  for  low-flux  imaging  scenarios  (as  in,  for  example,  observation  of 
dim  objects  in  space).  Here  we  must  again  rely  on  atmosphere-invariant  observables  like  the 
bispectrum  and  closure  phase.  Many  approaches  have  been  developed  for  estimation  of  the 
Fourier  phase  from  the  closure  phase  relations.  Several  early  papers  on  the  subject  (Rogstad, 
1968),(Bartelt  et  al,  1984)  proposed  solving  these  equations  recursively.  Namely,  two  initial 
Fourier  phases  6  a  and  0g  associated  with  two  connected  baselines  A  and  B  can  be  initialized 
arbitrarily  without  loss  of  generality  (see  Section  3.5.1).  The  Fourier  phase  9c  of  the  residual 
baseline  C  closing  A  and  B  can  be  determined  as:  9q  =  Pabc  ~  where  Pabc  is 

the  closure  phase  associated  with  triangle  ABC.  These  initial  phase  estimates  can  then 
bootstrap  solution  for  the  remaining  phases  using  the  appropriate  closure  relations  involving 
redundant  copies  of  the  baselines  associated  with  these  initial  phases.  While  simple  to 
implement,  such  approaches  are  clearly  sub-optimal  in  terms  of  estimation  accuracy;  the 
fact  that  information  about  a  given  Fourier  phase  is  contained  in  multiple  closure  phases 
implies  there  is  a  benefit  in  joint  inference  of  the  Fourier  phases  from  the  complete  set  of 
closure  phases. 

As  we  will  quantify  later  in  this  Section,  closure  phases  are  corrupted  with  a  variable 
amount  of  measurement  noise;  closure  phases  involving  short-baseline  measurements  will 
typically  incur  less  distortion  due  their  high  visibility  (and  hence  high-SNR)  relative  to  their 
long-baseline  counterparts.  Reliable  joint  inference  from  the  closure  phases  then  entails 
minimization  of  a  weighted  sum  of  closure  or  bispectrum  residuals  (i.e.  weighted  least- 
squares  methods).  Such  methods  are  of  course  complicated  by  wrapping  issues  analogous  to 
those  described  in  Chapter  3.  To  address  these,  Haniff  (1991)  developed  an  algorithm  which 
respects  the  27r-periodicity  of  the  closure  phase.  This  approach  uses  a  conjugate-gradient 
routine  to  find  the  set  of  Fourier  phases  {0}  which  minimizes  the  following  non-linear 
objective: 


Yi  = 


Nc 

E 


i= 1 


(  mod  2n{fii,cl  ~  (fl;i  +  @i2  + 


2 


(4.2) 
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where  Nc  is  the  number  of  closures,  p,-/C/  is  the  z’-th  closure  phase,  the  On,  0l2/  and  d&  are  the 
Fourier  phases  associated  with  this  closure,  and  Wi  is  a  weighting  factor  proportional  to  the 
estimated  variance  of  the  closure  measurement. 

While  this  minimization  has  been  shown  to  be  successful  in  certain  cases,  it  is  prone  to 
stagnation  at  local  minima  due  to  the  non-smoothness  of  the  mod  (x)  function  (Negrete- 
Regagnon,  1996),(Thiebaut  and  Giovannelli,  2010).  An  alternative  least-squares  approach 
suggested  by  Gorham  et  al.  (1989)  which  instead  works  with  the  closure  phasor  (i.e.  the 
normalized  bispectrum)  has  hence  generally  been  preferred  (Negrete-Regagnon,  1996).  This 
approach  replaces  the  objective  above  with  one  taking  the  form: 

2 

N  _  g/(0;i+0j2+0,3) 

T*  =  E - ^ -  <43> 

Objective  Y2  is  clearly  a  continuous  function  and  hence  local  minima  can  be  reliably  found 
with  gradient-descent  techniques.  However  the  objective  is  non-convex  in  general,  and  hence 
gradient-based  techniques  are  not  guaranteed  to  converge  to  the  true  global  minimum. 

To  the  best  our  knowledge,  work  by  Lannes  and  Anterrieu  (1999)  was  the  first  to 
rigorously  study  the  use  of  closure  phases  within  the  RSC  framework.  In  this  work,  the 
authors  present  a  mixed,  integer  linear-least  squares  formulation  of  the  phase  estimation 
problem  in  which  the  integer  component  pertains  to  the  unknown  (integer)  phase  wraps. 
They  first  propose  the  use  of  techniques  from  lattice  theory  to  resolve  the  wrap  ambiguities. 
They  then  develop  a  pseudo-inverse  estimator  which  recovers  the  Fourier  phases  from 
the  measured  baseline  phases  via  an  intermediate  computation  of  the  closure  phases. 
The  computation  of  the  necessary  covariance  matrices  for  optimal  estimation  within  this 
formulation  is  omitted.  Since  it  is  in  fact  the  closure  phases  which  are  the  direct  observables 
in  practical  low-flux  scenarios,  we  present  instead  an  approach  working  solely  with  the 
closure  phases  and  compute  the  necessary  covariances  in  straightforward  manner.  Lannes 
(2003)  later  proposed  an  alternate  estimator  based  on  computation  of  the  Smith  Normal  Form 
of  the  matrix  mapping  Fourier  phases  to  closure  phases.  This  estimator  has  the  advantage 
of  completely  eliminating  the  effect  of  phase  wrapping  in  the  closure  measurements  when 
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appropriate  routines  from  lattice  theory  are  used  to  pre-process  the  phase  measurements. 
In  contrast,  our  estimator  reliably  uses  standard,  fast  linear  estimation  techniques,  thereby 
obviating  computation  of  the  Smith  Normal  Form  at  the  possible  expense  of  an  extra  shift 
in  the  recovered  image.  Such  image  shifts  are  anticipated  to  be  of  negligible  importance  in 
most  cases  of  practical  interest. 

Though  the  bispectrum  and  associated  closure  phase  have  been  standard  interferometric 
observables  for  decades,  it  is  not  difficult  to  imagine  situations  in  which  use  of  these 
observables  unnecessarily  limits  reconstruction  performance.  In  Figure  4.1,  we  illustrate 
this  idea  with  an  example.  We  label  each  baseline  of  an  interferometric  array  with  a 
visibility,  which  is  an  indicator  of  the  strength  of  its  Fourier  component  relative  to  the 
overall  brightness  of  the  image.  In  accordance  with  the  power-law  decay  of  intensity  with 
spatial-frequency  modulus  in  an  overwhelmingly-large  fraction  of  natural  images  (see,  e.g. 
Ruderman  (1994)),  we  consider  a  visibility  distribution  which  drops  sharply  with  baseline 
length.  In  standard  bi-spectral  imaging,  the  high-spatial-frequency  phase  information 
associated  with  long  baseline  b*  is  recovered  via  forming  closures,  for  example,  with  a  short 
baseline  along  with  another  long  baseline  ( Traditional ).  As  we  will  see  later  in  the  chapter,  for 
a  certain  range  of  photon  fluxes,  the  SNR  of  the  bispectrum  is  approximately  proportional 
to  the  sum  of  the  reciprocals  of  the  squared  visibilities  of  the  associated  baselines  (Kulkarni 
et  al,  1991),  i.e.: 


(3  N 

l-2 

«=1  it  , 


(4.4) 


where  n  is  the  number  of  photoelectrons  (pe)  received  per  interference  fringe  per  exposure 
frame.  As  we  will  show  in  this  chapter,  this  SNR  model  extends  in  the  natural  way  to 
the  SNR  of  higher-order  observables  like  that  shown  on  the  right  ( Generalized ).  In  contrast 
with  the  Traditional  observable.  Generalized  observable  utilizes  only  short  (high-visibility) 
baselines  to  close  the  long  baseline,  resulting  in  a  near-doubling  of  the  resulting  SNR.  We 
will  henceforth  refer  to  the  higher-order  generalization  of  the  bispectrum  triple  product  as 
the  n-spectrum,  and  its  phase  as  the  generalized  closure  phase. 
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Figure  4.1:  Generalizing  the  phase  closure  concept .  SNRs  given  assume  each  aperture  contributes  h  —  2e3 
photoelectrons  to  each  fringe  in  pairwise  combination. 


Generalization  of  the  closure  phase  is  not  a  novel  concept  in  the  astronomical  community. 
Recently  Martinache  (2010)  proposed  observables  known  as  Kernel  phases,  which  are  obtained 
by  computing  an  orthogonal  basis  for  the  left-nullspace  (e.g.  using  the  Singular  Value 
Decomposition)  of  the  matrix  A  in  Equation  (4.1)  .  The  vectors  of  this  basis  are  assembled 
as  the  rows  of  a  matrix  K.  When  this  matrix  is  applied  to  the  measurements  f,  the 
contributions  due  to  the  atmosphere  (i.e.  <p)  clearly  cancel,  leaving  only  linear  combinations 
of  the  object  phases  in  9.  These  linear  combinations  can  then  be  used  to  reconstruct  the 
Fourier  phase  via  direct  application  of  the  pseudo-inverse  of  K  (Martinache,  2014),  or  via 
regularized  inversion  using  prior  image  models  (Ireland,  2013).  These  studies  consider 
masked-aperture  or  segmented-pupil  scenarios  in  which  either  closure  or  Kernel  phases 
can  be  obtained  directly  via  linear  combinations  of  the  raw  phase  measurements.  Since 
we  consider  the  scenario  encountered  in  long-baseline  interferometry,  in  which  accurate 
raw  phase  measurements  are  not  reliably  obtained  via  multi-frame  averaging,  our  inference 
problem  must  be  formulated  in  terms  of  integrated  phasor  observables,  i.e.  the  n-spectra 
defined  above. 
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Our  approach  is  also  conceptually  similar  to  the  approach  of  baseline  bootstrapping 
introduced  in  Chapter  3,  Section  3.5.3.  To  see  this,  let  us  define  the  normalized  optical 
path  difference  (OPD)  as  the  difference  in  the  atmosphere  phases  between  a  given  pair  of 
apertures.  In  the  bootstrapping  approach  first  proposed  by  Roddier  (1988),  temporal  OPD 
changes  along  a  long,  low-SNR  baseline  b  are  inferred  by  combining  (or  bootstrapping)  those 
measured  on  several  short  high-SNR  baselines  along  a  closed  path  involving  b.  In  principle, 
the  phases  on  all  baselines  can  then  be  de-rotated  accordingly  and  the  corresponding 
phasors  integrated  coherently,  thereby  boosting  the  SNR  on  all  baselines.  This  tracking 
and  de-rotation  can  be  done  real-time  with  fringe  tracking  or  offline  on  processed  data 
(Buscher,  2015).  As  in  bootstrapping,  our  approach  seeks  to  build  reliable  proxy  observables 
for  the  phase  estimates  of  the  low-visibility  baselines  by  forming  closed  paths  involving 
high-visibility  baselines.  However,  bootstrapping  eliminates  differential  phase  relative  to 
a  reference  exposure  frame  and  hence  implicitly  requires  a  threshold  per-frame  SNR  on 
the  shortest  baselines  for  robustness.  In  contrast,  our  approach  has  no  such  threshold  as  it 
works  directly  with  atmosphere-invariant  observables. 

In  this  chapter,  we  develop  a  comprehensive  algorithmic  framework  for  image  recon¬ 
struction  based  on  RSC  techniques  and  sparse  recovery,  and  analyze  its  performance  for 
the  two  popular  beam-combining  architectures  used  in  optical  interferometry  ( pairwise ,  and 
Fizeau).  The  chapter  is  organized  as  follows.  In  Section  II,  we  derive  an  SNR  model  for  the 
n-spectrum,  as  well  as  for  the  covariance  amongst  distinct  n-spectra.  The  analysis  begins 
with  the  pairwise  case  for  simplicity,  and  mathematical  arguments  for  an  extension  to  the 
Fizeau  case  are  presented  at  the  end  of  the  section  and  in  the  Appendix.  In  Section  III,  we 
describe  a  systematic,  integer-least-squares  approach  for  reconstructing  imagery  from  object 
visibilities  and  generalized  closure  phases  (GC's).  We  draw  on  the  notion  of  a  minimum 
cycle  basis  to  select  a  minimum  set  of  linearly-independent  GC's  (i.e.  a  basis )  of  minimum 
variance  which  spans  the  subspace  of  all  closures.  Leveraging  techniques  first  suggested  by 
Lannes  and  Anterrieu  (1999)  and  the  notion  of  wrap-invariant  measurement  mappings,  we 
apply  algorithms  from  lattice  theory  to  reliably  unwrap  the  GC's.  We  then  quantify  the  gain 
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in  theoretical  performance  afforded  by  generalizing  the  classical  third-order  atmosphere- 
invariant  observables  to  higher-order  cycles,  and  corroborate  this  theoretical  performance 
with  simulation  results.  In  this  simulation,  we  apply  our  algorithm  to  reconstruct  a  dim, 
structured  object  from  generalized  closures.  For  the  final  image  reconstruction,  we  employ  a 
sparse  recovery  algorithm  which  enforces  smoothness  in  the  image  domain  by  minimizing 
a  total  variation  metric. 

4.3  Preliminaries 

In  this  chapter  we  analyze  the  performance  of  generalized-closure-based  imaging  for  both 
the  pairwise  and  Fizeau  interferometric  architectures  (see  Section  1.2)  in  terms  of  mean- 
squared  error  (MSE)  of  the  phase  estimates.  We  begin  with  the  pairwise  architecture  since 
it  is  significantly  easier  to  analyze.  We  then  present  a  series  of  approximations  which 
accurately  describe  the  algorithm's  performance  in  the  Fizeau  case  while  keeping  the 
tediousness  of  the  mathematics  at  a  manageable  level.  In  both  cases,  the  analysis  is  based 
on  the  computation  of  the  moments  of  the  Poisson  distribution;  in  the  pairwise  case,  the 
first  two  moments  are  required,  whereas  in  the  Fizeau  case,  the  first  o  moments  are  needed, 
where  o  is  the  order  of  the  generalized  closure. 

4.3.1  Fringe  Noise  Model  for  Pairwise  Beam  Combiner 

Given  Nap  apertures,  recall  that  the  Fourier  phase  and  amplitude  of  the  object  is  encoded  in 
a  series  of  ( )  interference  fringes  observed  on  a  focal  plane.  In  the  pairwise  case,  each 
fringe  is  measured  on  a  separate  set  of  detectors  so  that  the  Poisson  shot  noise  incurred 
is  independent  for  each  fringe.  Suppose  each  aperture  receives  n  photons  and  that  this 
light  is  split  evenly  (e.g.  by  beam-splitter)  before  being  combined  with  the  light  from  the 
other  apertures.  Hence  each  aperture  contributes  n  =  N  "  ,  photons  to  each  fringe.  The 
expectation  of  number  of  photons  q  at  a  given  detector  k  for  a  given  fringe  is  given  by: 
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(4.5) 


<?(*))  =  ^(1  +  7  cos{cok  +  9  +  (p )) 

where  Np  is  the  number  of  pixels  over  which  these  photons  are  spread,  7  is  the  visibility  of 
the  fringe  (which  takes  it  value  between  0  and  1),  cv  is  the  fringe  frequency,  9  is  the  Fourier 
phase,  and  tp  is  the  atmosphere  phase. 

To  retrieve  the  object  phase  and  amplitude  from  its  fringe  encoding,  we  simply  take  the 
Fourier  transform  and  evaluate  it  at  the  fringe  frequency: 

1 

z  =  £  q(k)e~icok  (4.6) 

k= 0 

The  expectation  of  this  quantity  is  given  by: 

(z)  =  (4.7) 


where  is  the  phase  of  the  fringe. 

To  obtain  the  bispectrum  SNR,  it  will  be  useful  to  compute  the  power  in  the  fringe,  i.e. 
expectation  of  the  quantity  zz* : 


<22*>  =  (  E  E 

\  k=0  k'—Q 


This  is  equivalent  to: 


(zz*)  =  E  E  <#)?(*')> 

k= 0  k’=0 


g—iwkgiwk’ 


Using  the  first  two  moments  of  the  Poisson  distribution,  we  can  write: 


(4.8) 


(4.9) 


(q(k)q(k'))  =  (q(k)){q(k'))+5w{q(k)) 


(4.10) 


Substituting  into  Equation  (4.9),  we  obtain: 


Np— 1 


Np  — 1 


Np— 1 


(ZZ*)  =  £  <#)>*-“■*  £  w +  £  m) 
k= 0  V= 0  k= 0 


Simplifying  we  obtain: 


(zz*)  =  7  2n2  +  2n 


(4.11) 


(4.12) 
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4.3.2  N-Spectrum  Covariance  and  SNR  Models  for  the  Pairwise  Architecture 

Let  us  define  a  generalized  N -spectrum  G  as  the  product  of  complex  baseline  measurements 
along  an  o-edge  cycle.  Given  variances  crl(G)  and  ajm  (G)  of  real  and  imaginary  components, 
respectively,  we  define  the  pseudo-variance  of  the  product  as: 


Vvairunse(G)  :  =  (T2Re  +  ojm  =  (GG*)  -  (G)(G*) 


(4.13) 


or: 


VpairwiseyG)  — 


n  m) 


i=  1 


-n  <*.■><*?> 


(4.14) 


i= 1 


^pairwise  (  G  )  — 


FI (»27  i  + 2 n) 

i= 1 


-«2orK 


(4.15) 


Z=1 


In  the  specific  case  of  the  classic  triple  product  (i.e.  the  bispectrum),  this  expands  to 
(Kulkarni  et  al,  1991): 


Vpairwise(G)  =  +  7273  +  737i)  +  4n4(jf  +  72  +  73)  +  8” 


,.2„,2 


2„.2\ 


(4.16) 


The  SNR  Spairwise  :=  Y/^  is  then  given  by: 


5 


717273H2 


pairwise 


(4.17) 


\/  >!2(7?72  +  72 72  +  727i)  +  2/1(7?  +  7§  +  72)  +  4 

To  gain  insight  into  the  limiting  behavior  of  the  SNR,  we  make  two  approximations: 
A  low-SNR  approximation  can  be  found  by  considering  only  the  constant  term  in  the 
denominator,  which  will  dominate  when  n  and/ or  the  7,-  are  small. 


5 


pairwise, low 


717273^ 


(4.18) 


A  high-SNR  approximation  can  be  found  by  considering  the  highest-order  term  in  the 
denominator,  which  will  dominate  when  n  is  large  and/or  the  7 ,■  are  large. 


5 


pairwise, high 


7i 7273«2 


»2(7?722  +  7227|  +  7|7?) 


(4.19) 
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By  averaging  the  fringe  phasors  over  a  sufficient  number  of  frames  Nf,  we  can  build 
the  SNR  to  a  level  at  which  the  bispectrum  phase  (i.e.  the  closure)  phase  can  be  measured. 
For  phasor  SNRs  sufficiently  greater  than  1,  we  can  apply  the  following  well-known 
approximation  for  the  variance  on  the  corresponding  phase  9cj  (see  e.g.  Walkup  and 
Goodman  (1973)): 


~  NfS2 

We  then  obtain  the  following  approximations  for  the  closure  phase  variances 
SNR,  we  have: 

2  _ ,  4 

ed,low  Ny7j7273«3 

Taking  the  logarithm  we  find  that  the  variance  decouples: 

log  aedJozv  ~  -2  log  71-2  log  72-2  log  73-3  log  n  +  C  (4.22) 

where  C  is  a  constant  independent  of  the  flux  level  n  and  the  visibilities  7 
At  high  SNR,  we  have: 

ali’hish  ~  Wf  ^  +  +  ^ )  (4'23) 

Generalizing  the  results  above  to  n-spectra,  we  obtain  for  low  SNR, 


(4.20) 
.  At  low 

(4.21) 


GL./oiu 


Nfn°  YIU  7 


(4.24) 


and  for  high  SNR, 


%,,high  ~ 

where  0  is  the  order  of  the  generalized  closure. 
Equation  (4.26),  we  obtain: 


(4.25) 


Note  that  if  we  take  the  logarithm  of 


2  0 

l°g  gL low  ~  0  !°g  V  -  2  £  log  7 i  -  log  Nf 
n  i=  1 


(4.26) 
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4.3.3  Covariance  Matrix  of  the  Generalized  Closure  Phases  for  the  Pairwise 
Architecture 

In  this  section  we  generalize  results  of  Kulkarni  et  al.  (1991)  to  obtain  expressions  for 
the  n-spectra  covariance  matrices.  We  will  first  compute  the  covariance  between  two  n- 
spectrum  phasors  Gi  =  Re[Gi]  +  jIm[Gi\,  with  E[Gi]  =  gi,  and  G2  =  Re[G2]  +jIm[G2], 
and  E[G2]  =  g2-  Let  I  and  /  denote  the  set  of  baselines  in  n-spectra  G\  and  G2,  respectively. 
Then  we  have: 


(r{Gi,G*2) 


iei  jej  /  iei  jej 


(4.27) 


We  can  then  factor  out  the  factors  common  to  G\  and  G2  from  those  which  are  distinct 


to  obtain: 


cr(G!,G|) 


n  ( ^4)-  n  (^{4)  n  &)  n  (4) 

Jcein/  keinj  j  iei\J  mej\i 


Making  the  appropriate  substitutions  from  Equation  (4.12),  we  obtain: 


(4.28) 


n  (n27fc  +  2n4u,z2,J  -  H  n2j2k 

keinj  keinj 

x  n  ”7/  n  ”7m  (4-29) 
lei\J  mej\i 

where  W,  k~2k  =  1  if  the  common  measurements  are  the  same,  and  zero  if  they  are  conjugates 
of  each  other.  Similarly, 


cr(Gi,G|)  =  e2^Gi-ZG2) 


E[  («27fc  +  2n4u,z2,J-  H  fi2y2k 
keinj  keinj 

x  II  ”7/  n  ”7m  (4-30) 

lei\J  mej\i 

where  S2]  k:Zjk  =  1  if  the  common  measurements  are  conjugates  of  each  other,  and  zero  if 


a{G1,G1)  =  eln>^+^ 
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they  are  the  same. 

We  can  now  compute  real  and  imaginary  n-spectra  covariance  as: 


cr(Re[Gi],  Re[G2\)  =  ^Re[Cov(Gx,G^)  +  Cov{G1,G2)}  (4.31) 

cr(hn[Gi]rRe[G2])  =  ^lm[Cov{G\,  G|)  +  Cov(Gi,G2)}  (4.32) 

a(Re[Gi],  hn[G2\)  =  ^Im[— Cov(Gi,  G|)  +  Cou(Gi,  G2)]  (4.33) 

cr(Re[Gi\,  Im[G2\)  =  ^Re[-Cov(Glr  G|)  -  Cov(G1,G2)]  (4.34) 


Given  the  n-spectra  covariance  expression  above,  we  can  now  approximate  the  covariance 
between  two  generalized  closure  phases  using  Taylor  series  approximation  of  the  moments, 
which  is  a  technique  known  as  the  delta  method  (see,  e.g.  Casella  and  Berger  (2002)).  First  we 
define  the  angle  function  of  a  phasor  as  9{x,y)  =  arctan  The  corresponding  generalized 
closure  phases  are  given  by:  ©i  =  9(Im[Gx],Re[Gi})  and  ©2  =  9(hn[G2],Re[G2])-  Applying 
the  delta  method  yields  the  following  approximation: 


(4.35) 

+  (|l8l)  (|jJ^MGi],R<-[G2])  (4.36) 

+  (||J  (|lJ"('mlGll',m[G2|)  (437) 

(438) 

Evaluation  of  the  derivatives  in  the  expression  above  yields: 
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cr(©i,02) 


-Im\gl 


Re[g2 


+ 

+ 


Mg iP  +  Imigi}2J  VRe[£2p  +  lm[gi]2 

Re[gl\  \  (  - Im[g2 } 


Re[g i]2  +  Im[gl]2J  V^fe]2  +  Im[g2]2 
Mgi]  \  (  Mgi] 

Re[g i]2  +  Im[gi]2J  VRe[g2]2  +  Im[g 2}2 
~Im\gi\  \  (  ~Im[g2 


+ 


Mgi}2  +  Im[gi]2J  \Re[g2}2  +  Im[g2 } 


cr(Re[Gi\,  Im[G2]) 
a(Im[Gi\,  Re[G2\) 
a(Im[Gi\,  Im[G2\) 

cr(Re[Gi],  Re[G2\)  (4.39) 


4.3.4  Covariance  Matrix  Approximations  for  the  Fizeau  Architecture 

In  the  Fizeau  case.  Equation  (4.12)  becomes: 

(zz*)  =  7  2n2  +  Napii  (4.40) 

As  shown  in  Appendix  C.l,  we  can  approximate  the  variance  of  a  n-spectrum  of  order  o 


as: 


b Fizeau  (G) 

Hence  the  SNR  can  be  approximated  as: 


n  + n«ph) 

i= 1 


-«2‘Tb? 

i= 1 


(4.41) 


5, 


Fizeau 


VimnU'Yi 


(4.42) 


[n?=i  (n^f  +  Napn)]  -n^YlLni 
Applying  the  same  approximations  in  Appendix  C.l  to  the  covariance  yields  the  follow¬ 
ing  analogues  to  Equations  (4.29)-(4.30) 


cr(Gi,G|)  =  eZTT/^Cn-ZGz) 


x 


n  (»27fc  +  KpnSZ]  krZ2rk)  -  n  n2y2k 


Ikeinj 


keinj 


n  nrYi  Ft  n7m  (4-43) 

lei\J  mej\l 
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and. 


a{GlrG2)  =  e^i^G1+ZG2) 


x 


n  (^Tic  +  NapnSZ]  krZ2ik)  -  n  »27 k 
keinj  keinj 


x  n  ”7/  n  n7m 

/6l\/  mej\l 


(4.44) 


4.4  Fourier  Phase  Recovery  using  the  N-Spectra 

In  this  section,  we  briefly  review  the  integer  least-squares  approach  for  RSC  phase  recovery 
(see  e.g.  Lannes  and  Anterrieu  (1999),  Kurien  et  al.  (2016)  for  more  comprehensive  treat¬ 
ments),  and  then  describe  an  algorithm  for  selection  of  a  set  of  generalized  closure  relations 
of  minimum  total  variance. 


4.4.1  RSC  Phase  Recovery  with  Wrap-Invariant  Closure  Mappings 

We  begin  by  recalling  the  fundamental  RSC  equation  from  Chapter  3: 


M 


9 

$ 


=  (5  +  2rce 


(4.45) 


In  this  Section  we  discuss  generalized  closure  relations  which  eliminate  the  (p  contribution 
in  the  model  in  Equation  (4.45).  We  first  recall  from  Chapter  3  (see  Definition  3.5.1)  that 
following  the  analysis  of  Lannes  and  Anterrieu  (1999),  the  range  of  the  matrix  M  can  be 
divided  into  two  subspaces:  (1)  the  subspace  K  spanned  by  the  d  columns  associated  with 
the  Fourier  phases,  or  equivalently  the  range  of  the  sub-matrix  Mg  formed  by  these  columns, 
and  (2)  the  piston,  or  aberration,  phase  space  L  is  the  span  of  the  Nlip  columns  associated 
with  the  piston  phases,  or  equivalently  the  range  of  the  sub-matrix  M,^  formed  by  these 
columns.  Note  that  the  d  columns  spanning  K  are  linearly-independent  by  virtue  of  having 
non-zero  entries  in  mutually  disjoint  sets  of  baseline  indices.  Hence  dim(K)  =  d.  The 
subspace  L  has  dimension  Nap  —  1  (Lannes  and  Anterrieu,  1999);  the  constant  vector  forms 
the  one-dimensional  nullspace  of 
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Additionally,  let  us  establish  the  following  Definitions: 

Definition  4.4.1.  (Interferometric  graph):  The  interferometric  graph  of  an  array  is  the 
directed  graph  whose  vertices  are  the  Nap  apertures  in  the  array  and  edges  are  the  (^p) 
baselines  connecting  all  vertex  pairs. 

Definition  4.4.2.  (Cycle-space  of  a  directed  graph)  (Liebchen  and  Rizzi,  2005):  Given  a 
directed  graph  G  with  a  set  of  vertices  V  and  edges  E,  the  cycle  space  of  the  graph  is 
the  vector  space  of  QIeI  spanned  by  the  incidence  vectors  of  cycles  with  a  (clockwise,  or 
counter-clockwise)  orientation.  It  can  be  shown  that  the  dimension  of  the  cycle  space  is 
M  +  N  —  1,  where  M  =  |  E  |  and  N  =  \  V  \ . 

Suppose  we  have  a  basis  for  the  cycle-space  of  the  interferometric  graph.  Let  us  stack  the 
elements  of  this  basis  as  row  vectors  in  a  (' Nap2  a)  x  (n|p)  matrix  Cmc,  which  maps  baseline 
phase  measurements  to  generalized  closure  phases.  Recall  that  closure  relations  eliminate 
piston  differences  in  the  measurements  so  that  Cmc  annihilates  the  subspace  L,  i.e.  the  space 
spanned  by  the  columns  of  M  corresponding  to  <p. 

In  Figure  4.2,  we  show  a  simple  four-aperture  interferometric  graph  and  one  possible 
cycle  basis  for  the  graph.  For  this  particular  example,  we  can  write: 


C 


me 


11  1  0  0  0 
00-1110 
10  0  011 


(4.46) 


where  the  column  (i.e.  baseline)  indexing  in  the  matrix  follows  the  labeling  in  the  Figure. 
Note  that  the  final  triangle  cycle  (i.e.  the  one  containing  baselines  2,  4,  and  6,  and  represented 
by  the  row  vector  w*  :=  (0, 1,0, 1,0,  —1)),  can  be  represented  as  a  linear  combination  of  the 
cycles  in  the  basis:  w*  =  ci  +  C2  —  C3,  where  {c,}  are  the  rows  of  Cmc 

While  for  purposes  of  simplicity  we  have  chosen  a  cycle  basis  consisting  exclusively  of 
three-baseline  cycles,  in  general  the  elements  of  a  cycle  basis  can  contain  any  number  of 
edges. 


Lemma  4.4.1.  The  nullspace  of  Cmc  is  L. 
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Figure  4.2:  One  possible  cycle-basis  for  a  simple  interferometric  graph.  Note  the  fourth  triangle  (i.e.  baseline 
set  {2,4,6})  can  be  expressed  as  a  linear  combination  of  the  cycles  shown. 


Proof :  Clearly  the  Nnp  columns  of  L  are  in  the  nullspace  of  Cmc  from  the  arguments  above. 
To  see  that  these  columns  span  the  nullspace,  note  that  from  the  well-known  Rank-nullity 
theorem  from  linear  algebra  (see,  e.g.  Bretscher  (2001)),  we  have: 


dim(ker( Cmc))  =  ^  =  Nap  -  1  (4.47) 

Hence  dim(ker( Cmc))  =  dim{L).  □ 

We  can  now  form  closure  relations  as: 


6 

0 

= 

C  oc  0 

$ 

<P 

(4.48) 


where  Coc  :=  CmcM q,  and  is  the  mapping  between  object  phases  and  closure  phases. 


Proposition  4.4.2.  Let  Arx  and  Ar^  denote  the  vectors  containing  the  x-  and  y -coordinates  of  the 
baselines  in  an  array ,  respectively.  If  the  array  is  a  valid  RSC  array ,  these  columns  form  a  basis  for 
the  two-dimensional  nullspace  of  Coc. 


Corollary  4.4.3.  For  a  valid  RSC  array,  the  mapping  Coc  is  injective  up  to  an  image  shift. 


Proofs  of  Proposition  4.4.2  and  Corollary  4.4.3  are  given  in  Appendix  C.2.  □ 

Let  us  define  yc;  as  the  observation  vector  of  wrapped  generalized  closure  phases,  and 
2necj  the  wrapping  vector  (with  integer  ec/).  Our  phase  measurement  model  is  then  given 
by: 
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Coc9  =  yP  +  n 


(4.49) 


where  ye  =  y ci  +  27rec/,  and  n  is  the  measurement  noise,  which  we  shall  assume  in  this 
chapter  is  predominantly  due  to  shot  noise  on  the  focal  plane. 

Given  that  the  generalized-closure  mapping  is  injective  up  to  an  image  shift,  we  can 
now  formulate  the  well-posed  recovery  of  the  Fourier  phases  as  the  following  generalized, 
integer  least-squares  problem,  which  finds  a  vector  in  the  range  of  the  closure  matrix  at  a 
minimum  Mahalanobis  distance  from  the  actual,  wrapped  measurement  vector: 

Orsc  =  ar§mined/e{ye  -  Coed)  £_1  (ye  -  C0,f)  (4.50) 

where  E  is  an  estimator  of  the  covariance  matrix  of  the  phase  measurements.  This  estimator 
can  be  obtained  by  substituting  estimates  for  the  object  visibilities  {7/}  into  the  covariance 
expressions  in  the  previous  section.  It  is  acknowledged  that  this  least-squares  solution  is 
neither  the  optimal  nor  maximum-likelihood  solution  to  the  phase  inference  problem.  For 
one,  a  Gaussian  distribution  on  the  closure  phases  is  implicitly  assumed  to  hold,  and  can 
only  approximately  hold  for  high-SNR  closure  phases.  It  is  therefore  assumed  that  n-spectra 
will  be  integrated  for  a  sufficiently  large  number  of  frames  for  this  approximation  to  be 
reliable  for  most,  if  not  all,  of  the  closure  phases. 

Assuming  the  covariance  matrix  estimator  is  not  degenerate,  it  will  admit  a  Cholesky- 
decomposition  E  =  BBT.  Equation  (4.50)  is  then  equivalent  to  searching  for  vectors  e  which 
minimize  the  projection  of  a  whitened  measurement  B  1  ye  onto  the  subspace  defined  by 
ker((B~1Coc)T).  Specifically  we  seek  to  minimize: 

/(e)  =  ||PzB-1(yc/+27rec/)||2  (4.51) 

where  Py  is  a  matrix  representing  the  orthogonal  projection  from  1R"  onto  ker((B  ]Coc)T). 

Letting  e'd  =  —eci,  we  can  rewrite  the  above  objective  function  as: 

/(e')  =  ||PzB-1(yc/-27re';)||2  =  ||PEB-1yd-27rPsB“1e';||2  (4.52) 
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This  optimization  problem  is  equivalent  to  the  so-called  closest  vector  problem  (CVP)  in 
the  theory  of  lattices.  The  connection  of  interferometric  phase-unwrapping  with  the  CVP 
problem  was,  to  the  best  of  our  knowledge,  first  explored  in  Lannes  and  Anterrieu  (1999).  In 
this  paper,  the  authors  formulate  an  analogous  minimization  problem  for  the  raw  baseline 
measurements. 

We  will  define  a  lattice  L(Z")  as  the  set  of  points  generated  by  integer  combinations 
of  the  column  vectors  of  a  matrix  L.  Letting  P  =  PvB  \  our  optimization  problem  then 
amounts  to  the  following:  Find  the  lattice  point  in  P(Z")  which  is  closest  to  Pyc;/  (2tt).  A 
compact  representation  of  the  lattice  T  is  given  by: 

{m<n—(d+N— 3)  1 

flW  |  Vflf  ezl  (4.53) 

where  { v, }  are  linearly-independent  and  together  form  a  basis  of  the  lattice.  Given  the 
lattice  basis,  several  algorithms  exist  for  finding  the  closest  lattice  point  to  a  specified  vector. 
A  popular  class  of  algorithms,  known  as  the  Sphere-Decoding  algorithms,  are  efficient 
searches  for  the  closest  lattice  point  within  a  hypersphere  of  a  certain  radius  centered  on  the 
input  vector  (see  e.g.  Agrell  et  al.  (2002)).  For  the  simulations  in  this  chapter,  we  instead  use 
the  lower-complexity  Babai  Nearest  Plane  (Babai-NP)  algorithm  (Babai,  1986).  For  lattice 
bases  which  are  nearly  orthogonal  (such  as  those  we  use  for  our  simulations),  this  algorithm 
offers  reliable,  albeit  not  guaranteed,  performance  in  practice. 

Suppose  we  have  found  a  basis  for  the  lattice  P(Z”),  and  we  have  solved  the  Closest 
Vector  Problem  for  a  given  closure  vector  yc/.  Let  b  be  the  output  of  the  Babai  Nearest 
Plane  Algorithm  -  i.e.  it  is  the  lattice  point  which  is  the  closest  to  y^.1  We  now  seek  to  solve 
for  the  wrap  vector  corresponding  to  this  lattice  point,  i.e.  we  seek  a  solution  to: 

b*  =  Pe  (4.54) 

ht  is  well-known  that  the  performance  of  the  Babai  algorithm  and  other  CVP  algorithms  is  improved  when 
the  lattice  basis  is  reduced  as  mentioned  in  Section  3.3.1.  We  used  the  LLL  algorithm  implementation  due  to 
Zhou  (2014)  in  our  simulations. 
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Note  that  P  is  a  projection  matrix  and  thus  not  full-rank,  and  therefore  there  will  be 
infinitely-many  solutions  to  this  equation.  Indeed,  by  analogy  to  the  integer  wrap  vectors 
in  Chapter  3  (which  were  correct  only  up  to  an  integer  vector  in  the  range  of  M),  our 
solution  yd  is  correct  only  to  within  an  integer  vector  eresict  in  the  range  of  Coc.  Hence  there 
is  a  fundamental  ambiguity  that  needs  to  be  addressed.  Once  the  measurement  vector  is 
unwrapped,  there  are  various  ways  to  solve  the  RSC  system  which  respect  the  presence  of  a 
residual  ambiguity  eresiy  described  above.  Following  standard  least-squares  principles,  we 
first  compute  via  projection  the  vector  y*d  uniC}  in  im( Coc)  closest  (in  Mahalanobis  distance) 
to  the  unwrapped  measurement  vector,  i.e. 


yeW>  =  B"1(I-p£)By'i  <4-55> 

Per  the  first  method  (Lannes,  2003),  we  can  compute  the  Smith  Normal  Form  (SNF)  (see 
Theorems  3.5.2  and  3.5.3)  {U,  D,V}  of  our  matrix  Coc,  and  set  Uc  =  U-1,  Dc  =  D,  and 
Vc  =  V  1 .  We  then  can  write: 


Coc  =  UcDcVc  (4-56) 

A  valid  RSC  solution  can  then  be  obtained  by  evaluating: 

eRSC  =  Vc'D+U  cVcUmlQ  (457) 

where  Djt  denotes  the  pseudo-inverse  of  Dc- 

The  effect  of  the  residual  ambiguity  vector  eres^  on  this  solution  requires  careful  consid¬ 
eration.  Note  that  the  resulting  wrap-induced  error  in  this  solution  is  given  by: 

2neRSC  =  27rVc1D+Ucj1eres!rf  (4.58) 

Since  the  matrices  Uc  and  Vc  are  unimodular,  they  are  invertible  over  the  integers.  If  all 
the  elementary  divisors  in  Dc  are  equal  to  1,  then  will  also  be  integral,  and  hence  eRsc 
will  be  integral.  This  in  turn  guarantees  that  the  final  error  in  the  Fourier  phase  will  be  0 
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mod  2n,  i.e.  that  the  RSC  solution  9r$c  is  wrap-invariant. 

Definition  4.4.3.  A  wrap-invariant  closure  mapping  Coc  is  one  whose  elementary  divisors 
of  the  matrix  Coc  are  all  1. 

Proposition  4.4.4.  (Sufficient  condition  for  wrap-invariant  closure  imaging):  Given  a  wrap- 
invariant  closure  mapping  Coc,  the  RSC  phase  solutions  derived  from  the  associated  (generalized) 
closure  phases  will  he  immune  to  the  effects  of  phase-wrapping. 

Recall  that  in  Chapter  3,  we  observed  that  this  condition  can  be  violated  by  RSC  patterns 
which  are  not  wrap-invariant.  Conversely,  we  have  observed  that  with  wrap-invariant 
patterns,  the  closure  mappings  associated  with  the  cycle  bases  that  we  will  consider  in 
this  paper  (i.e.  fundamental  cycle  bases  and  minimum  cycle  bases)  are  also  typically 
wrap-invariant  in  practice. 

As  it  turns  out,  given  a  wrap-invariant  closure  mapping,  the  inverse  problem  in  Equation 
(4.49)  can  be  solved  reliably  even  without  the  Smith  Normal  Form;  it  is  sufficient  to  find  any 
unimodular  r-by-r  sub-matrix  C  of  Coc  and  solve  the  associated  smaller  system  in  Equation 
(4.59)  to  obtain  a  valid  solution.  Recall  that  for  valid  RSC  arrays  Coc  will  have  r  +  2  columns 
and  hence  in  this  approach,  two  (non-collinear)  object  phases  are  implicitly  set  to  zero.  This 
implicit  selection  then  fixes  the  fundamentally-ambiguous  translation  of  the  scene  discussed 
at  the  beginning  of  the  section. 


0  =  C-yclMC)  (4.59) 

Since  C  is  unimodular,  C  1  will  have  solely  integral  entries  so  that  the  resulting  wrap 
error  C  1  eresui  =  0  mod  2tt.  Assuming  the  measured  phase  has  been  correctly  unwrapped, 
the  solution  in  Equation  (4.59)  amounts  to  a  so-called  basic  solution  of  our  generalized 
least-squares  problem.  We  then  have  the  following  expression  for  the  error  covariance 
matrix  for  the  estimator  9rsc  (K ay,  1993). 

Zbnsic  =  (CrE“1C)“1  (4.60) 
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The  basic  solutions  belong  to  the  countably-infinite  set  of  solutions  to  Equation  (4.49).  It 
turns  out,  somewhat  surprisingly,  that  any  particular  solution  within  this  set  can  be  reliably 
limited  to  a  mere  image  shift  for  patterns  that  are  wrap-invariant.  The  complete  set  is  given 
by: 


=  C+y*clMC)  +  0O  (4.61) 

where  do  is  any  vector  in  the  nullspace  of  Coc,  and  C+  denotes  the  pseudo-inverse  of  Coc. 
The  pseudo-inverse  can  be  computed  using  the  Singular- Value-Decomposition  (SVD)  of  Coc, 
which  is  given  by: 


Coc  =  UaZaVa  (4.62) 

where  Ut7  and  Vo-  are  (Nflp2  x)  x  (N^  a)  and  d  x  d  orthogonal  matrices,  respectively.  X,r  is  a 
(N"p-i)  x  ^  diagonal  matrix  with  r  non-zero  diagonal  entries  (the  so-called  singular  values  of 
Coc  )/  where  r  =  rank( Coc)  =  d  —  2. 

The  Moore-Penrose  pseudo-inverse  is  then  given  by  (Bretscher,  2001): 

C+  =  UoE+Vj  (4.63) 

where  E(|  is  a  diagonal  matrix  whose  r  non-zero  diagonal  entries  are  the  reciprocals  of  the 
corresponding  non-zero  entries  in  Zp. 

In  Chapter  3,  we  showed  that  the  effect  of  wrapping  on  the  RSC  pseudo-inverse  solution 
to  Equation  (4.45)  can  be  reliably  limited  to  an  image  shift  for  wrap-invariant  patterns.  An 
analogous  result  and  proof  apply  to  closure-based  RSC  imaging: 

Proposition  4.4.5.  For  wrap-invariant  closure  mappings,  the  error  induced  by  wrapping  of  the 
(generalized)  closure  phases  can  be  limited  to  an  image  shift  using  the  standard  Moore-Penrose 
pseudo-inverse  estimator. 

Proof:  See  Appendix  C.3  □ 

We  can  leverage  this  result  to  establish  an  analogous  result  pertaining  to  the  well-known 
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non-linear  least-squares  formulation  due  to  Gorham  et  al.  (1989)  discussed  in  Section  4.2. 
Namely  we  have  the  following  Corollary: 

Corollary  4.4.6.  Given  a  wrap-invariant  closure  mapping,  the  set  of  solutions  minimizing  Objective 
Y2  in  Equation  4.3  differ  from  the  true  solution  by  an  image  shift  in  the  noiseless  case. 

Proof:  For  Objective  Y2  to  be  zero,  clearly  each  term  (i.e.  each  squared-residual)  in  the 
summation  must  be  zero.  Hence  we  must  have  e^''d  =  e/Yn  +^a+tfa)  |  y j  These  constraints 
are  clearly  equivalent  to  the  linear  phase  constraints  solved  by  the  family  of  solutions  in 
Equation  (4.61).  Hence  the  solution  sets  are  the  same  in  the  noiseless  case.  □ 

The  error  covariance  of  elements  in  the  pseudo-inverse  solution  can  be  expressed  as 
(Montgomery  et  al.,  2006): 


Zpnw  =  VjZ2+W  (4.64) 

where  Vcr  is  obtained  by  omitting  the  final  2  columns  of  V(r  which  form  a  basis  for  the 
nullspace  of  Coc. 

4.4.2  Selection  of  the  N-Spectra 

In  this  section,  we  describe  a  strategy  for  obtaining  a  near-optimal  linearly-independent  set 
of  generalized  closure  phases.  This  strategy  is  founded  on  the  notion  of  minimum  cycle  basis 
from  graph  theory. 

Definition  4.4.4.  ( Minimum  cycle  basis  of  a  directed  graph):  Let  each  edge  of  the  graph  be 
assigned  a  positive  weight  cy,  and  the  weight  of  a  cycle  be  defined  as  the  sum  of  the  weights 
of  its  constituent  cycles. 

With  these  Definitions,  we  can  then  ask  for  the  set  of  linearly-independent  oriented 
cycles  which  spans  the  cycle  space  and  has  minimum  total  weight.  We  call  such  a  set  the 
minimum  cycle  basis  of  a  directed  graph  G. 

We  next  state  a  fundamental  lemma  that  establishes  a  practical  method  for  obtaining  a 
minimum  cycle  basis,  for  which  the  proof  can  be  found  in  many  textbooks  on  graph  theory 
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(see  e.g.  Gross  and  Yellen  (2006)). 


Lemma  4.4.7.  (Optimality  of  the  Greedy  Algorithm):  Consider  a  vector  space  V,  a  set  of  vectors 
which  span  V,  and  a  corresponding  set  of  weights  for  each  of  these  vectors.  A  basis  of  minimum  total 
iveight  can  be  found  by  first  sorting  the  vectors  in  decreasing  order  of  weight,  and  then  selecting 
vectors  for  the  basis  in  this  order  if  and  only  if  they  are  linearly-independent  of  previous  selections. 
This  is  known  as  the  greedy  algorithm. 

From  Lemma  4.4.7,  we  can  find  a  set  of  linearly-independent  closure  relations  of  minimum 
total  variance  by  employing  the  greedy  algorithm  on  the  set  S  of  all  possible  closure  relations. 
However,  this  set  is  equivalent  to  the  set  of  all  possible  cyclic  permutations,  which  becomes 
enormous  for  arrays  of  large  size.  Fortunately,  in  an  extension  of  a  result  due  to  Horton 
(1987)  for  undirected  graphs,  Liebchen  and  Rizzi  (2005)  showed  that  the  elements  of  the 
minimum  cycle  basis  of  a  directed  graph  could  be  found  among  a  special  (and  much  smaller) 
subset  of  S  known  as  the  Horton  cycles.  We  describe  these  cycles  through  the  following 
definitions: 

Definition  4.4.5.  ( Shortest  path  tree):  The  shortest  path  tree  of  a  graph  G  with  respect  to  a 
node  A  is  the  spanning  tree  which  connects  node  A  to  each  other  node  in  G  via  a  path  of 
minimum  weight. 

Figure  4.3  gives  a  simple  example  of  a  shortest  path  tree  for  an  interferometric  graph  in 
which  the  top-most  node  is  the  root.  The  construction  of  shortest  path  trees  is  a  well-studied 
problem  in  graph  theory  (see  e.g.  Gross  and  Yellen  (2006))  which  is  typically  solved  using  a 
shortest-path  routine  such  as  Dijkstra's  Algorithm  (Dijkstra,  1959). 

Definition  4.4.6.  (Horton  cycle):  Given  the  shortest  path  tree  T  originating  from  an  arbitrary 
node  A,  a  Horton  cycle  is  a  cycle  formed  by  connecting  the  endpoints  of  any  two  branches 
of  T. 

Lemma  4.4.8.  (Minimum-cycle-basis  elements  are  Horton  cycles):  All  elements  of  the  minimum- 
cycle-basis  are  Horton  cycles.  Therefore  in  searching  for  the  minimum-cycle-basis,  it  suffices  to  search 
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0.3 


-  -  -  Shortest-Path  Tree 


Figure  4.3:  A  shortest  path  tree 

the  Horton  cycles.  Proof:  We  summarize  the  Proof  due  to  Liebchen  and  Rizzi  (2005)  in  Appendix 

C.4.  □ 

To  exploit  the  Lemma  above,  we  seek  a  decomposition  of  the  weight  of  a  cycle  into  the 
weights  of  its  constituent  edges  (i.e.  baselines).  As  the  following  Propositions  show,  at 
the  extreme  low-  and  high-SNR  limits,  the  n-spectrum  variance  does  indeed  allow  such  a 
decoupling. 

Proposition  4.4.9.  (Loiv-SNR  Decoupling  Approximation):  In  the  lozv-SNR  limit,  the  logarithm 
of  the  total  variance  of  any  cycle  is  given  (up  to  a  global  additive  constant)  by  the  sum  of  the  costs 
assigned  to  each  baseline  in  the  cycle,  inhere  the  cost  of  a  baseline  k  with  visibility  7*.  is  given  by: 
ck  =  logf  -21og7fc. 

Proof:  This  can  be  seen  by  re-writing  the  low-SNR  approximation  in  Equation  (4.26)  as: 

log  al,,low  ~  E  lo§  \  -  2  lo§  7«  +  C  (4.65) 

i= 1  n 

where  C  =  log  Nj  is  a  constant  with  respect  to  both  visibility  and  flux  and  therefore 
irrelevant  for  the  purposes  of  closure  selection.  □ 

Proposition  4.4.10.  (High-SNR  Decoupling  Approximation):  In  the  high-SNR  limit,  the  total 
variance  of  any  cycle  is  given  (up  to  a  global  scaling  factor)  by  the  sum  of  the  costs  assigned  to  each 
baseline  in  the  cycle,  where  the  cost  of  a  baseline  i  with  visibility  7 ,  is  given  by:  c,  =  ^j. 

Proof:  This  follows  directly  from  the  high-SNR  approximation  in  Equation  (4.25).  □ 
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Based  on  the  Propositions  above,  we  now  have  a  way  to  construct  a  set  of  linearly- 
independent  generalized  closures  of  minimum  variance  at  the  low-  and  high-  SNR  extremes. 
As  demonstrated  in  the  next  section,  even  at  moderate  flux-levels  where  these  approxi¬ 
mations  lack  precision,  this  closure  selection  method  can  yield  significant  improvements 
over  traditional  closure  imaging  restricted  to  baseline  triangles.  To  estimate  the  scaling  of 
the  computational  cost  of  the  method,  recall  that  the  greedy  algorithm  is  run  on  the  set  of 
Horton  cycles.  Since  there  are  Nap  shortest-path  trees,  each  rooted  at  a  distinct  aperture,  and 
m  —  ( Nap  —  1)  cycles  in  each  of  these  trees,  the  total  number  of  Horton  cycles  is  0(mNap). 
Since  m  =  0(N%p),  we  see  that  the  algorithm  is  0(N%p),  which  is  the  same  complexity  as  if 
the  greedy  algorithm  were  run  on  an  exhaustive  listing  of  all  three-baseline  closures. 


4.5  A  Practical  Algorithm  for  RSC  Closure  Imaging 

In  this  section  we  synthesize  the  techniques  developed  in  previous  sections  into  an  algorithm 
for  RSC  imaging  using  generalized  closures.  Recall  that  closure  selection  relies  on  estimates 
of  the  baseline  visibilities  {7;}.  These  estimates  can  be  derived  for  the  pairwise  case  by 
solving  Equation  (4.12)  for  7  to  yield: 


1/  (zz*)  —  2n 


(4.66) 


Similarly,  for  the  Fizeau  case,  we  can  solve  Equation  (4.40)  to  obtain: 


yj {zz*)  -  Napn 


(4.67) 


Having  estimated  both  Fourier  phases  and  visibilities,  we  can  form  estimates  for  the 
complex  visibilities  v  of  the  object.  To  estimate  the  image  coefficients  x  from  these  complex 
visibilities,  we  must  apply  a  deconvolution  algorithm.  For  this  purpose  we  have  chosen  the 
standard  sparse-recovery  regularization  known  as  Total  Variation  Minimization  (TV-Min), 
whose  success  in  reconstructing  piecewise  smooth  objects  is  well-known(Becker  et  ah,  2011). 
The  TV-Min  regularization  is  given  by: 
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(4.68) 


x  =  argmina  ||a||Ty  subject  to:  ||Fa  —  v||2  <  e 

where  F  is  a  partial  Fourier  matrix  whose  rows  are  vectorized  representations  of  the  2D 
sinusoids  associated  with  the  array's  measured  spatial  frequencies. 

Algorithm  3  provides  the  complete  set  of  steps  we  have  discussed. 


Algorithm  3  RSC  Image  Reconstruction  Algorithm 

1.  Estimate  Visibilities  7 ,  (c.f.  Equations  (4.66)-(4.67)) 

2  Select  N-Spectra  Relations  via  Minimum  Cycle  Basis 

2.1  choose  decoupling  approximation  according  to  SNR  regime  (c.f.  Proposition  4.4.9  or 
4.4.10) 

2.2  assign  baseline  edge  weights  accordingly 

2.3  enumerate  the  Horton  cycles  of  the  interferometric  graph 

2.4  execute  greedy  algorithm  on  set  of  Horton  cycles  to  find  minimum-variance  set  of 
m  —  ( Nap  —  1)  n-spectra 

3.  Average  the  Selected  N-Spectra 

4.  Unwrap  the  generalized  closure  phases  corresponding  to  these  N-Spectra  (c.f.  Sec¬ 
tion  4.4.1) 

5.  Solve  generalized  least-squares  problem  with  unwrapped  generalized  closures 

6.  De-convolve  the  estimated  complex  visibilities  to  produce  image  reconstruction 
using  the  regularization  in  Equation  (4.68)  or  other  regularization  technique 


We  used  the  software  package  known  as  NESTA  (Becker  et  al,  2011)  to  execute  regular¬ 
ization  in  Equation  (4.68). 


4.6  Algorithm  Performance 

In  this  Section,  we  present  the  results  of  application  of  Algorithm  3  to  a  simulated  scenario 
of  imaging  a  structured  object  space. 

4.6.1  Sensitivity  Limits 

A  useful  benchmark  for  our  results  is  provided  by  the  Cramer-Rao  Lower  Bound  (CRLB) 
for  interferometric  phase  estimation  in  the  absence  of  atmospheric  turbulence.  We  begin  by 
defining  our  fringe  measurement  model  in  the  standard  way: 
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p  =  A 


y  cv 


nN"P-nT 

N%  N}r 


(4.69) 


where  p  is  the  vector  of  pixel  counts,  xcz„  and  ycv  are  the  (N“p)  real  and  imaginary  components 
of  the  complex-visibility  components,  respectively,  A  is  the  matrix  mapping  the  real  and 
imaginary  parts  of  the  complex  visibilities  to  pixels  on  the  focal  plane  (often  called  the 
visibility-to-pixel  matrix  V2PM  in  the  literature),  and  1^2  is  the  all-ones  vector  of  length 
N2p.  Recall  that  in  Section  1.4  we  formed  this  mapping  as  a  matrix  A,  whose  rows  were 
the  sinusoidal  functions  associated  with  each  fringe  generated  by  the  beam  combiner.  A 
comparison  between  the  definitions  of  the  two  mappings  yields  the  following  relation: 


A  = 


In 

w. 


AP, 


(4.70) 


where  Pa  performs  a  permutation  of  the  columns  of  A  to  respect  the  ordering  of  the 
quadrature  components  xcv,  and  ycv  in  Equation  (4.69). 

Leveraging  the  analysis  in  Chapter  1  (see  Section  1.4)  2,  we  can  compute  the  Fisher- 
Information  Matrix  (FIM)  for  complex-visibility  estimation  with  a  Fizeau  beam-combiner 
as: 


Wc  =  (A^A)  (4.71) 

where  is  a  diagonal  matrix  whose  diagonal  entries  are  the  expected  values  of  the  photon 
counts  at  each  detector  in  the  array. 

The  CRLB  for  the  covariance  of  these  vector  parameters  is  then  given  by: 

cw„  -  >  o  (4.72) 

where  the  notation  >  in  this  context  means  that  the  matrix  difference  on  the  left-hand  side 
is  positive-semidefinite. 

2The  interested  reader  is  directed  to  Zmuidzinas  (2003)  for  the  original  derivation 
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The  above  bound  applies  to  single-frame  fringe  phasor  estimation  for  a  non-redundant 
array.  Since  we  are  instead  evaluating  phase  estimation  schemes  for  a  redundant  array 
given  multiple  frames  of  data,  we  must  incorporate  these  transformations  into  the  bound 
(Kay,  1993).  Namely,  for  a  redundant  array,  the  parameter  vector  must  be  shortened  to  the 
d  distinct  complex-visibility  phasors  xcvj  +  (lj)ycv,d-  Let  us  define  the  vector  zcv^  as  the 
concatenation  of  xCV/d  and  y cv,d-  The  fringe  model  is  then  given  by:  p  =  ARzczv/  +  , 

where  R  is  a  2(N"P)  x  2d  matrix  mapping  the  real  and  imaginary  parts  of  the  distinct  complex 
visibilities  to  those  of  the  2(N"f ’)  generated  fringes.  The  FIM  then  becomes: 


i*„,  =  (RtatA7;ar 


(4.73) 


Defining  the  function  gg  =  arctan  the  CRLB  for  the  phase  covariance  matrix  "Lg  can 


be  expressed  as: 


1  d§6  (zcv,d)  j— 1  dgg(z Cv,d)  q 

'  dzr  - 


Nframes  ^Z 


cv,d 


'■‘cv.d 


(4.74) 


where  is  the  Jacobian  matrix  of  the  distinct  Fourier  phases  with  respect  to  the 


distinct  complex-visibility  phasors. 

Since  the  left-hand-side  of  Inequality  (4.74)  is  positive  semi-definite  (PSD),  and  the 
diagonal  entries  of  a  PSD  matrix  must  be  non-negative,  we  arrive  at  the  following  useful 
bound  for  the  variance  of  the  estimated  phases: 


var(6j)  > 


1 


N 


frames 


dge(z  CV,  rf)rl  dgg (z  cv,d  ) 
d%cv,d  Ll,d  dzcv(j 


(4.75) 


4.6.2  Simulation 

In  this  section,  we  provide  the  results  of  a  simulation  in  which  Algorithm  3  was  applied 
to  both  interferometric  architectures.  For  these  simulations  we  used  an  RSC  pattern  of  the 
Y-pattern  type  as  shown  in  Figure  4.4.  The  corresponding  UV-sampling  for  this  pattern 
is  displayed  in  Figure  4.5.  To  prevent  aliasing  on  the  simulated  focal  plane,  this  aperture 
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pattern  was  mapped  onto  the  Golay  non-redundant  pattern  (Golay,  1970)  shown  in  Figure 
4.6  for  fringe  generation;  the  model  emulates  the  redundant-to-non-redundant  pupil  re¬ 
mapping  technique  developed  by  Perrin  et  al.  (2006)  and  described  in  the  Introduction. 
Our  simulation  assumed  Poisson-distributed  shot-noise  and  idealized  detectors  with  zero 
read  noise,  and  so  is  representative  of  a  shot-noise-dominated  scenario.  For  bispectrum 
observables,  we  used  the  unbiased  estimator  reported  by  Gordon,  J.  A.  and  Buscher,  D.  F. 
(2012): 


glib  =  ZaZbZc  -  \za\2  ~  \zb\2  -  |zc|2  +  2 Nnpn  (4.76) 

where  za,  zb,  and  zc  are  the  fringe  phasors  associated  with  the  three  sides  of  a  bispectrum 
triangle,  and  n  is  an  estimate  of  the  total  number  of  photons  per  aperture  incident  upon  the 
focal  plane.  3 

The  target  selected  was  NASA's  Cloud- Aerosol  Lidar  and  Infrared  Pathfinder  Satellite 
Observations  (CALIPSO)  satellite  for  which  the  truth  image  (Hill,  2008)  is  shown  in  Figure 
4.7.  Figure  4.8  shows  the  image  at  the  resolution  attainable  by  the  pattern.  Two  flux 
levels  (2000  photoelectrons /aperture /frame,  and  500  photoelectrons /aperture /frame)  were 
considered  for  these  simulations.  Comparative  error  analysis  using  Equation  (4.60)  showed 
a  lower  predicted  Root-mean-squared  (RMS)  phase  error  for  the  decoupling  approximation 
in  Proposition  4.4.10  than  that  in  Proposition  4.4.9,  and  hence  the  former  was  chosen  in 
Algorithm  3.  Bispectra  and  n-spectra  were  integrated  for  5e4  frames,  which  corresponds 
approximately  to  an  8-minute  observation  time  if  we  assume  a  typical  frame  duration  of  10 
msec.  An  implementation  of  the  Bellman-Ford-Moore  shortest-path  algorithm  (O'Connor, 
2012)  was  used  to  generate  a  minimum  cycle  basis  as  per  Algorithm  3. 

To  show  the  potential  impact  of  generalizing  the  closure  phase  notion.  Algorithm  3  was 
compared  with  an  analogous  algorithm  which  instead  used  the  minimum-variance  set  of 
(Np  1  j  ^dependent  traditional  (i.e.  three-baseline)  closure  phases.  This  minimum  set  was 

3We  did  not  derive  the  more  complicated  bias  corrections  for  higher-order  nspectra  as  empirical  results 
indicated  that  at  the  flux  levels  considered,  the  importance  of  these  bias  corrections  declined  with  the  order  of 
the  nspectra. 
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Figure  4.4:  RSC  aperture  pattern  used  in  simidation 
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Figure  4.5:  UV -sampling  for  RSC  pattern 
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Figure  4.6:  Go/czi/  non-redundant  beam-combiner  pattern 
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Figure  4.7:  Truth  image  for  simulation:  the  CAL1PSO  satellite 
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Figure  4.8:  Truth  image  at  the  resolution  of  the  interferometric  pattern 

generated  via  application  of  the  greedy  algorithm  described  above  to  the  set  of  all  ( 'Nap2f 1). 

Numerical  results  for  the  pairwise  and  Fizeau  architectures  in  the  higher  flux  scenario 
are  shown  in  Figures  4.9  and  4.10,  respectively.  In  these  plots,  the  RMS  Fourier-phase 
errors  for  10  independent  simulation  trials  are  plotted  as  a  function  of  their  corresponding 
visibilities  for  the  basic  solution.  Median-filtered  versions  of  the  predicted  standard  errors 
based  on  Equations  (4.60)  and  (4.64)  are  also  shown  for  both  basic  and  pseudo-inverse 
solutions,  respectively.  4 

As  expected,  the  Fizeau  architecture  proves  to  be  more  sensitive  than  the  Pairwise 
architecture,  which  corroborates  analysis  by  Zmuidzinas  (2003).  Figures  4.11  and  4.12 
respectively  show  Fizeau  performance  for  the  higher  flux  scenario  at  a  shorter  integration 
time  (lc3  frames),  and  for  the  lower  flux  scenario  at  the  long  integration  time  (5c4  frames). 
Note  that  in  the  use  of  independent  generalized  closure  phases  according  to  Algorithm 
3  outperforms  the  analogous  scheme  in  which  only  traditional  closures  are  used;  the 

4The  actual  predictions  exhibit  oscillations,  which  reflect  the  particular  structure  of  our  mapping  between 
Fourier  phases  and  GC's.  For  clarity  we  have  created  smoothed  versions  of  the  predicted  curves  by  creating 
bins  with  edges  at  visibilities  {5, 10, 20, 40, 80, 160, 320}  *  10~3,  and  displaying  the  median  values  in  each  bin. 
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root-mean-squared  (RMS)  errors  are  noticeably  lower  for  the  former. 

It  is  interesting  to  examine  performance  relative  to  the  CRLB  described  in  Section 
1.4.  Recall  that  our  CRLB  is  an  atmosphere-oracle  bound  and  hence  it  should  bound  phase 
estimation  accuracy  for  the  case  of  a  stationary  fringe  averaged  over  many  frames.  The 
variance  of  such  an  estimator  can  be  easily  calculated  by  leveraging  the  analysis  in  Section 
4.3.  Namely  the  pseudovariance  of  the  phasor  is  given  by  applying  Equation  (4.40)  after 
averaging,  so  that  SNR  of  the  averaged  fringe  phasor  z  is  given  by: 

cr?  =  1  =  _  (4  77) 

z  SNR 2  |j£[z]||2  2N framesn'Y2  V  ' 

We  compute  the  predicted  phase  RMS  for  the  atmosphere-oracle  by  taking  the  square- 
root  of  this  expression,  and  add  it  to  the  plots.  From  the  plots,  we  see  that  CRLB  is  virtually 
identical  to  this  expression,  except  for  isolated  dips  in  the  CRLB.  These  dips  correspond  to 
the  redundant  baselines,  which  as  expected,  admit  increased  phase  sensitivity  due  to  their 
multiplicity. 

In  the  higher-flux  scenario  (n  =  2e3),  our  algorithm's  phase-estimation  performance  is 
within  a  factor  of  2  from  the  CRLB  for  the  vast  majority  of  sampled  visibilities.  As  the  flux 
drops  to  n  =  5e2,  the  performance  begins  to  diverge  from  the  CRLB  as  the  closure  phase 
variance  leaves  the  regime  described  by  Equation  (4.25)  and  enters  the  regime  described  by 
Equation  (4.24).  That  is,  the  variance  is  no  longer  accurately  modeled  simply  by  a  linear 
combination  of  the  individual  phase  variances  of  the  form  in  Equation  (4.77);  rather  it 
becomes  inversely-proportional  to  the  product  of  the  squared-visibilities.  The  rapid  decrease 
in  fidelity  of  the  closure  phase  at  low-SNR  is  well-known  in  interferometry  (Kulkarni  et  al., 
1991). 

Sample  image  reconstructions  are  shown  in  Figure  4.13.  The  reconstructions  derived 
from  generalized  closures  show  greater  fidelity  to  the  true  image  than  those  derived  from 
traditional  closures,  and  thereby  corroborate  the  numerical  analysis  presented  in  the  plots. 

It  is  noteworthy  that  while  our  algorithm  attempts  to  find  a  minimum-variance  set  of 
n-spectra,  this  set  is  not  necessarily  the  optimum  set  for  phase  estimation;  our  algorithm 
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Figure  4.9:  Pairwise  Phase  Recovery  Results  for  Flux  n  —  2e3  pe/ap, /frame,  5e4  frames 
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Figure  4.10:  Fizeau  Phase  Recovery  Results  for  Flux  n  —  2e3  pe/ap/frame,  5e4  frames 
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Figure  4.11:  Fizeau  Phase  Recovery  Results  for  Flux  n  —  2e3  pe/ap/fmme,  le3  frames 
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Figure  4.12:  Fizeau  Phase  Recovery  Results  for  Flux  n  —  5e2  pe/ap/frame,  5e4  frames 
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Figure  4.13:  Fizean  Image  Reconstruction  Results.  The  reconstructions  in  the  left  column  used  traditional, 
three-baseline  observables,  whereas  those  in  the  right  column  used  generalized  closures  selected  according 
to  Algorithm  3.  (top  row)  n  —  2e3  photoelectrons/aperture/frame  (pe/ap/frame),  5e4  frames,  (middle  row) 
n  =  2e3  pe/ap/frame,  2e3  frames,  (bottom  row)  n  —  5e2  pe/ap/frame,  5e4  frames. 


107 


does  not  account  for  the  conditioning  of  the  matrix  Coc  resulting  from  a  particular  choice  of 
nspectra.  In  fact  there  may  be  some  cases  in  which  the  use  of  traditional  closures  results 
in  lower  phase-estimation  error  due  to  better  conditioning  in  the  corresponding  Coc.  In 
practice  the  decision  to  use  traditional  or  generalized  closures  can  be  resolved  by  examining 
the  phase  errors  predicted  by  Equation  (4.60),  for  basic  solutions,  or  Equation  (4.64)  for 
pseudo-inverse  solutions. 

4.6.3  Generalized  Closures  in  Non-Linear  Least  Squares  Approaches 

This  paper  has  proposed  both  novel  reconstruction  methodology  as  well  as  new  interfero¬ 
metric  observables.  A  natural  question  arises  as  to  whether  the  advantages  of  the  latter  are 
retained  with  other,  existing  methodology.  In  particular,  we  consider  the  class  of  techniques 
which  solve  the  non-linear  least  squares  inference  problem  involving  bispectrum  phasors 
(Gorham  et  al.r  1989),  (Negrete-Regagnon,  1996)  (see  Section  4.2). 

While  these  algorithms  bring  the  complexities  involved  in  non-linear  optimization  (e.g. 
possible  stagnation  in  local  minima  and/ or  slow  convergence),  they  also  feature  the  ability 
to  utilize  any  number  of  closure  relations.  Algorithm  3  restricts  attention  to  closure  relations 
which  form  cycle  bases,  thereby  guaranteeing  well-posed  Fourier  phase  recovery  while 
keeping  the  size  of  the  associated  CVP-unwrapping  problem  tractable.  However,  as  noted 
by  Kulkarni  et  al.  (1991),  fitting  to  all  ( )  closures  will  in  principle  improve  estimation 
performance  when  the  per-frame  flux  is  low,  due  to  the  fact  that  closure  phases  de-correlate 
as  the  flux  decreases  5.  Specifically,  as  computed  by  Kulkarni  et  al.  (1991)  for  traditional 
closures  in  the  pairwise  case,  the  correlation  coefficient  of  two  closures  sharing  a  common 
baseline  approaches  i  in  the  high-SNR  regime  (i.e.  for  fry2  3>  1)  as  one  would  expect  6.  On 
the  other  hand,  the  correlation  coefficient  in  the  low-SNR  regime  can  be  accurately  modeled 

1  4 

as  }i  «  |7  H  . 

To  assess  the  phase-estimation  performance  provided  by  generalized  closures  relative 

5The  interested  reader  is  directed  to  Buscher  (2015)  for  a  dedicated  discussion  of  this  issue 

6Here  Kulkarni  et  al.  (1991)  assumes  a  common  visibility  7  among  all  baselines  in  the  closure 
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to  that  using  the  complete  set  of  traditional  closures,  we  consider  the  algorithm  developed 
first  by  Gorham  et  al.  (1989)  which  works  the  normalized  bispectrum  instead  of  the  closure 
phase.  We  used  the  MATLAB®  non-linear  least-squares  solver  (NLS)  known  as  Isqnonlin 
(MATLAB,  2014)  to  minimize  the  objective  Y2  in  Equation  (4.3).  This  solver  obtains  a 
quadratic  approximation  of  the  objective  at  each  step  and  moves  towards  the  minimum  of 
this  approximation,  at  which  point  another  approximation  is  computed  and  the  process 
repeats.  The  weights  in  Y2  were  set  as  in  Gorham  et  al.  (1989): 


where  \\gi\\  is  the  magnitude  of  the  averaged  unbiased  nspectrum,  and  agi  is  the  empirical 
standard  deviation  of  its  quadrature  components. 

Results  of  applying  this  solver  to  the  lower-SNR  scenario  of  the  previous  section  are 
given  in  Figure  4.14.  A  single  randomly-chosen  initialization  point  was  used  in  all  cases. 
The  top  row  shows  reconstructions  and  associated  convergence  times  for  the  NLS  algorithm 
using  traditional  closures  (left)  and  generalized  closures  selected  from  a  minimum  cycle 
basis  (middle),  respectively.  The  bottom  row  gives  analogous  results  for  the  case  in  which 
all  (N^p)  traditional  closures  are  used.  Three  snapshots  at  iteration  counts  20,  40,  and  80  are 
shown.  Iterations  after  80  resulted  in  negligible  change  to  the  image  quality.  These  results 
suggest  that  generalized  closures  can  provide  at  least  the  sensitivity  of  traditional  closures. 
At  the  same  time,  they  benefit  from  the  increased  convergence  speed  afforded  by  a  problem 
size  proportional  to  the  size  of  cycle  basis  (i.e.  oc  N%p)  as  opposed  to  the  size  of  the  set  of  all 
closures  (i.e.  oc  N%p).  These  results  hence  show  that  generalized  closures  can  serve  as  more 
efficient  sources  of  phase  information  than  traditional  closures. 
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Convergence  in  5  sec  Convergence  in  4  sec  Convergence  in  3  sec 
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Figure  4.14:  (top  row)  Reconstructed  images  and  convergence  times  for  NLS  algorithm  with  traditional-closure  basis  (left),  generalized-closures  from 
minimum  cycle  basis  (middle),  and  generalized-closures  from  minimum-variance  spanning  tree  (right),  and  (bottom  row)  Reconstructed  images  and  elapsed 
running  times  with  all  (Nfp)  traditional  closures  at  iteration  20  (left),  iteration  40  (middle),  and  iteration  80  (right) . 


Figure  4.15:  Raw  ( left )  and  median-filtered  (right)  Fourier  phase  error  using  the  NLS  approach. 
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In  Figure  4.15,  we  plot  phase  error  as  a  function  of  visibility  for  the  NLS  approach  using 
a  cycle  basis  of  traditional  closures  (blue),  a  cycle  basis  of  generalized  closures  (green),  and 
using  all  (N^p)  traditional  closures  (red).  The  median  absolute  error  for  20  trials  is  shown  in 
the  scatter  plot  in  the  left  panel.  The  same  median-filter  that  was  applied  to  the  predictions 
in  Figures  4.9-4.12  was  then  applied  to  this  raw  data  to  produce  the  results  in  the  right  panel 
of  the  Figure.  It  is  clear  that  the  sensitivity  advantage  of  generalizing  the  closure  relations  is 
retained  in  the  NLS  approach  for  this  example. 

It  is  noteworthy  that  the  selection  of  a  minimum-variance  set  of  generalized  closures 
involves  computational  overhead  in  the  execution  of  the  greedy  algorithm.  Flowever,  it  is 
anticipated  that  in  many  cases  we  may  be  able  to  obtain  a  reliable  surrogate  for  this  set 
by  instead  by  forming  a  fundamental  cycle  basis  (Liebchen  and  Rizzi,  2005)  associated  with 
minimum-variance  spanning  tree  of  the  interferometric  graph.  Constructing  this  alternate 
basis,  which  was  introduced  in  Section  3.5.3,  obviates  the  need  for  execution  of  the  greedy 
algorithm;  it  merely  requires  the  execution  of  a  Minimum-weight  Spanning-tree  algorithm 
(e.g.  Prim-Jarnik  Algorithm  (Prim,  1957))  with  each  edge  weight  set  to  the  reciprocal  of  the 
squared  visibility  associated  with  baseline,  which  is  again  used  a  proxy  for  the  baseline's 
phase  variance  7.  The  elements  of  the  corresponding  minimum-variance  spanning  tree 

7We  used  a  MATLAB®  implementation  of  Prim's  algorithm  which  is  available  online  (Greenbaum,  2007) 
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Figure  4.16:  The  minimum-variance  spanning  tree 

(MVST)  cycle  basis  are  then  found  by  closing  the  edges  of  this  spanning  tree,  which  is  shown 
in  Figure  4.16.  Indeed,  many  of  the  Horton  cycles  comprising  the  minimum  cycle  basis 
computed  for  this  example  pattern  and  scenario  are  also  members  of  the  MVST  basis.  Not 
surprisingly,  the  reconstruction  with  this  basis  shown  in  the  upper  right  corner  of  Figure 
4.14  approaches  of  the  quality  that  for  the  minimum-cycle-basis  set  (compare  with  image  in 
the  middle  of  the  top  row  in  Figure  4.14).  Moreover  the  median  absolute  phase  error  for 
this  approach,  which  is  shown  in  cyan  in  the  right  panel  of  Figure  4.15,  is  very  similar  to 
that  observed  with  the  Minimum  Cycle  Basis. 

While  the  analysis  of  Kulkarni  et  al.  (1991)  revealed  the  advantage  of  supplementing  a 
cycle  basis  with  additional  traditional  closures  in  the  pairwise  case,  this  advantage  should 
in  principle  extend  to  the  case  of  generalized  closure  phases  in  the  Fizeau  architecture.  Our 
CRLB  comparison  above  (i.e.  Figures  4.10  and  4.11)  confirms  that  there  is  limited  scope  for 
improvement  in  the  regime  jj-  ~  60.  In  the  regime  explored  in  the  lower-SNR  scenario  (see 
Figure  4.12),  however,  there  is  certainly  a  CRLB-gap  which  may  be  closed  with  appropriate 
selection  of  supplementary  generalized  closures.  We  leave  the  development  of  algorithms  to 
perform  such  a  selection  efficiently  for  future  work. 
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4.7  Conclusion 


In  this  chapter,  we  have  developed  a  novel  method  for  interferometric  imaging  which 
employs  RSC  techniques  and  a  generalized  notion  of  bispectrum  observable  (the  n-spectrum). 
We  have  established  that  phase  estimation  from  these  observables  can  be  well-posed  for  valid 
RSC  arrays.  We  have  also  provided  a  fast  algorithm  for  selection  of  a  minimum-variance 
set  of  these  observables,  which  is  based  on  the  concept  of  minimum  cycle  basis  in  graph 
theory.  Then,  leveraging  the  lattice-theory  problem  formulation  first  proposed  by  Lannes 
and  Anterrieu  (1999)  for  unwrapping  of  the  closure  phases,  as  well  as  techniques  from 
sparse  recovery  for  image  reconstruction  in  the  presence  of  Fourier  undersampling,  we 
have  proposed  a  new  algorithm  for  image  reconstruction  in  optical  interferometry.  We 
have  shown  that  the  performance  of  the  phase-estimation  part  of  our  algorithm  can  be 
quantified  from  first  principles  using  standard  linear  estimation  theory.  We  have  used 
simulation  at  different  photon-flux  levels  to  corroborate  this  analysis,  and  to  show  the 
potential  advantage  of  performing  inference  using  the  n-spectrum  observable  as  opposed  to 
its  classical  counterpart  -  the  bi-spectrum.  It  is  our  hope  that  both  the  theoretical  framework 
employed  in  this  chapter,  as  well  as  the  practical  algorithm  itself,  will  prove  useful  to  future 
users  in  designing  RSC-based  interferometric  systems  and  processing  their  data. 
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Chapter  5 


Conclusions 


This  thesis  has  presented  both  theoretical  and  practical  contributions  to  the  study  of  robust 
imaging  in  optical  interferometry.  Our  principal  theoretical  contribution  is  that  we  have 
developed  the  notion  of  the  wrap-invariant  RSC  pattern.  We  have  shown  that  these  patterns 
enable  unique  recovery  of  the  Fourier  phase  of  the  object  under  observation.  These  results 
are  fundamental  in  the  sense  that  they  apply  to  the  entire  family  of  RSC  techniques  which  fit 
the  Fourier  phase  either  directly  to  the  complex  visibilities  measured  by  the  interferometer, 
or  atmosphere-invariant  combinations  of  these  such  as  the  bispectrum.  We  summarize  our 
uniqueness  results  in  Table  5.1,  classifying  them  according  to  the  Approach  type  (i.e.  Phase 
vs.  Phasor),  and  the  observable  used  in  the  phase  recovery  (i.e.  direct  measurements  of  the 
fringe  phase  or  phasor  vs.  an  atmosphere-invariant  combination  thereof).  References  to 
other  work  in  the  literature  employing  each  observable-approach  pair  are  also  provided  for 
context. 

We  have  leveraged  these  theoretical  results  to  construct  a  novel  algorithmic  framework  for 
image  reconstruction.  In  our  framework,  we  have  generalized  classical  atmosphere-invariant 
observables  to  improve  estimation  of  the  Fourier  phase  of  complex  objects.  The  ultimate 
goal  of  this  effort  was  to  achieve  near-immunity  to  the  effect  of  atmospheric  turbulence 
in  terms  of  Fourier-phase  estimation  accuracy,  or  more  concretely,  to  show  performance 
approaching  the  atmosphere-oracle  Cramer  Rao  Lower  Bound.  Using  our  algorithm,  we 
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provide  strong  theoretical  and  empirical  evidence  that  we  can  achieve  this  goal  in  the  regime 
in  which  the  per-frame  photon  flux  is  on  the  order  of  10  photons  per  interferometric  fringe 
(i.e.  jf-  is  on  the  order  of  10).  We  have  seen  that  the  increased  sensitivity  of  our  generalized 
observables  extends  to  other  existing  phase-recovery  algorithms  used  in  the  field.  Lastly,  we 
demonstrate  that  the  accurate  phase  estimation  provided  by  our  framework  allows  powerful 
techniques  from  Compressed  Sensing  to  generate  high-quality  interferometric  imagery. 
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Appendix  A 


Appendix  to  Chapter  1 

A.l  Sinusoidal  Dependence  of  Field  on  Interferometer  Focal  Plane 

As  an  aid  for  the  following  derivation.  Figure  A.l  shows  the  geometry  involved  in  beam 
combination  at  the  focal  plane. 

Note  that  we  can  express  the  vector  position  of  pixel  at  coordinate  p  relative  to  the 
position  of  aperture  j  as: 


Xj  =  d  +  p  -  r j  (A.l) 

Moreover  the  magnitude  of  the  wavevector  k ;  is  while  the  unit  vector  in  the  direction 
of  k j  is  given  by: 

*/  =  ,  ‘  „  ||2  ['/  +  0]  (A.2) 

Vlldll  +  llrdl 

Hence  the  phase  added  to  the  field  due  to  the  path  between  the  beam  combiner  and 
focal  plane  (at  pixel  p)  is  given  by: 

9  77 

ki '  xi  =  /  =  =  \rj  +  d]  •  [d  +  p  -  r j\  (A.3) 

AV  lldll  +llr;ll 

Expanding  the  dot-product  and  noting  that  d  •  p  =  0,  d  •  r;  =  0,  and  d  •  d  is  constant 
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3 


Figure  A.l:  Beam-Combination  Geometry 


independent  of  apertures,  we  obtain  (up  to  a  constant  offset  independent  of  aperture): 


k;  •  xi  = 


2n 


A 


I  j  II 2  ,  ||  1 1 2 

d  +  r,- 


lr; 


P)  +<Pj 


(A.4) 


where  <p;  := 


2n  r 


A^lldf+llpH2' 
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Appendix  B 


Appendix  to  Chapter  3 

B.l  Proof  of  Proposition  3.5.7 

In  this  Proposition,  we  propose  that  we  can  solve  for  the  Fourier  phase  vector  from  the 
closure  relations  in  two  separate  integer-preserving  steps:  the  first  ( Step  A)  involving  the 
SNF  decomposition  of  C m_>c,  and  the  second  (Step  B)  involving  that  of  M.  To  show  the 
validity  of  this  two-stage  decomposition,  we  require  the  following  Lemma. 

Lemma  B.1.1.  ker( COT_%c)  =  L 

Proof :  Clearly  every  vector  in  L  must  be  in  ker( Cm_>c),  since  otherwise  there  would  exist 
a  combination  of  {cp}  which  would  produce  a  non-zero  closure  phase.  This  is  impossible 
due  to  the  atmosphere-annihilating  property  of  the  closure  phase.  We  must  also  show 
that  L  spans  the  kernel.  Note  that  dim(im(C,n^c))  =  r  =  (N2  ')/  where  we  use  im(.)  to 
denote  the  range  of  a  matrix.  Noting  that  CmH>c  has  (^)  columns,  we  have  dim(ker (Cm_».c))  = 
Cl)  —  CC  ')  =  N  —  1  by  the  Rank-nullity  Theorem  from  linear  algebra.  Since  dim(L)  =  N  —  1, 
the  columnspace  L  indeed  spans  the  kernel.  □ 

Suppose  we  have  the  SNF  decomposition  of  CmH>c  =  UcDcVc-  The  result  of  Step  A  is 
then: 


h  =  Vc1Vcvc1(j  Cl¬ 


ine 


h,cl) 


(B.l) 
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Given  Lemma  B.1.1,  Equation  (B.l)  can  be  re-written  as: 


fa  =  fa  +  vl  +  27rV-'D+UcXd  (B.2) 

where  /3q  is  the  phase  measurement  corresponding  to  the  optimum  least-squares  solution 
for  Equation  (3.23),  and  vi  is  a  vector  in  L.  Since  the  elementary  divisors  of  Cm-,c  are  all 
1  by  construction  (by  Lemma  2.6),  the  rightmost  term  in  Equation  (B.2)  will  be  0  mod  2n 
and  hence  the  first  step  is  thus  integer-preserving.  Moreover,  since  V/  is  restricted  to  L,  the 
Fourier-phase  error  resulting  from  application  of  the  SNF-based  inverse  V^,1  D^U^,1  in  Step 
B  will  be  0  mod  2/r  if  all  elementary  divisors  of  M  are  1.  Hence  the  Proposition  holds.  □ 

B.2  Proof  of  Lemma  3.7.2 

In  this  section,  we  prove  Lemma  3.7.2,  which  states  the  following: 

Given  wrap-invariance,  the  (column)  vector  UxM^e;*  has  integer  entries  belozv  row  2. 

Given  that  we  have  a  wrap-invariant  pattern,  we  know  that  the  elementary  divisors  of 
M  are  all  1.  Hence  there  exists  an  integer  vector  ko  such  that: 

4  =  Mk0  (B.3) 

Substituting  Equations  (C.45)  and  (3.26)  into  Equation  (C.38),  we  obtain: 

etr  =  V^L+UjU^Vjko  (B.4) 

Noting  that  LG  is  orthogonal,  this  equation  can  be  simplified  to 

e<r  =  (V,  -  N)Vjk0  (B.5) 

where  N  is  a  matrix  of  the  same  size  as  \a,  which  is  zero  except  for  the  last  three 
columns.  These  last  three  columns  are  identical  to  those  of  Va,  and  hence  by  Lemma  4.1 
comprise  an  orthogonal  basis  for  the  nullspace  of  M.  Noting  the  orthogonality  of  V,r,  this 
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can  be  further  simplified  to: 


e„  =  k0  -  NVjko  (B.6) 

To  proceed,  the  following  Definition  will  be  useful: 

Definition  A.l  (The  canonical  basis  for  the  nullspace  of  M):  The  canonical  basis 
{w/},  i  G  1,2,3  for  the  three-dimensional  nullspace  of  M  can  be  derived  trivially  from  the 
well-known  tilt-position  degeneracy  in  interferometry  described  in  Section  3.5.1  (Wieringa, 
1992).  Namely,  we  can  define  the  basis  as  the  columns  of  a  (N  +  d)  x  3  matrix  W /a.f ( M :  as 
follows: 


Wfer( M)  =  {  W2  j  W2 1 W3}  = 


0 

Inxi 


Arx  Ary 


(B.7) 


where  rx  and  rv  are  the  x—  and  y— positional  coordinates  of  the  apertures  associated 
with  each  row,  respectively,  and  Arv  and  Arv  their  respective  pairwise  differences.  □ 

Note  that  each  of  the  three  non-zero  column  vectors  { v^},  k  G  1, 2, 3  in  N  can  be  expressed 
as  a  linear  combinations  of  the  elements  of  the  canonical  basis  {w ,},z  G  1,2,3  defined  above, 
i.e. 


Vfr  =  fljWi  +  fl2W2  +  «3W3  (B.8) 

As  in  Section  3.5.1,  let  us  again  use  K  to  denote  the  set  of  indices  in  e,r  associated  with 
the  Fourier  phases  (as  opposed  to  the  piston  phases),  and  their  corresponding  rows  in  N. 
Hence  we  have: 


Uxe(r,K  =  Uxk0/K  -  UxNxVjk0  (B.9) 

Since  Ux  is  an  integer  matrix,  the  first  term  is  clearly  integral.  Let  us  then  examine  the 
second  term,  and  in  particular,  the  product  UxN «.  By  substitution  from  Equation  (C.49),  we 
see  that  for  any  of  the  three  non-zero  columns  of  Nx,  we  have  (by  simply  changing  the 
basis  for  the  nullspace): 
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UxVjyc  =  UlUxWi^  +  fl2UxW2,K  +  03UXW3JC 


(B.10) 


where  { Wj  are  the  vectors  comprising  the  upper  partition  of  Equation  (B.7).  The  first 
term  in  Equation  (C.51)  is  trivially  0  since  is  the  zero-vector.  Now  note  from  Theorem 
3.5.2  that  Ux  is  a  matrix  which  annihilates  all  the  spatial  frequencies  in  the  matrix  X  below 
row  2.  But  from  Definition  A.l,  these  spatial  frequencies  are  identically  the  contents  of  the 
two  columns  and  (up  to  a  uniform  scaling  factor).  Therefore  the  column  vector  in 
Equation  (C.51)  is  zero  below  row  2.  This  means  in  turn  that  the  second  term  in  Equation 
(C.50)  is  zero  below  row  2,  and  hence  that  Uxe^x  is  integral  below  row  2  (since  the  first 
term  is  integral).  □ 
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Appendix  C 


Appendix  to  Chapter  4 


C.l  Fizeau  Variance  Approximations 

C.l.l  Variance  Decomposition 

In  this  section  we  justify  our  approximation  for  the  variance  of  the  n-spectrum  for  the  Fizeau 
architecture.  Our  goal  here  is  not  to  compute  all  of  the  numerous  terms  in  the  variance,  but 
rather  to  present  the  mathematical  intuition  behind  the  approximation  we  have  used.  The 
general  formula  for  the  variance  of  an  n-spectra  G  is: 

-< 

For  simplicity,  we  analyze  the  standard  bispectrum  case  ( o  =  3).  The  analysis  extends 
naturally  to  the  n-spectrum.  Note  that  by  extension  of  Equation  (4.8)  the  first  term  can  be 
written  as: 


(C.l) 


E  (Cl(pn)^(Pb)^(Pc)q(Pd)^Pe)q(Pf)) 

Pa,Pl„Pc,Pd,Pe,Pf 


X  eivi (pd-pa) eiu>2(Pe~Pb) eiV3(Pf~Pc)  (Q2) 
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Similarly,  the  second  term  in  Equation  (C.l)  can  be  written  as: 


(S*)  (fr* )  - 

Pa,Pb,Pc,Pd,Pe,Pf 


x  eivl(pd-pa)eiU2(Pe-Pb)eiU3(Pf-Pc)  (Q3) 


Proceeding  analogously  to  Equation  (4.10),  we  utilize  the  general  formula  for  the 
moments  of  a  Poisson  distribution  to  perform  the  following  decomposition  of  the  first 
factor  in  the  summand  in  Equation  (C.2)  (Kulkarni  et  al.,  1991): 
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(q(pa)q(pb)q(pc)q(pd)q(pe)q(pf))  =  (c.4) 

+{q(pa)){q(pb))(q(pc)){q(pd))(q(pe)){q(pf))  (c.5) 

+£p=p.=pAq(p))(q(Pb))(q(pc))(q(pe))(q(pf))  (c.6) 

+sP=ph=pe{q(p))(q(pa))  (q{pc))  (q(pd)){q(pf))  £7) 

+6p=prpf(q(p))(q(Pa))(q(Pb))(q(Pd))  (q(Pe))  (C.8) 

+6p=pa=Pd(q(p))sP'=pb=Pe(q(p'))(q(pc))(q(Pf ))  (C9) 

~^^p=pa=pd(q(p))^p'=pc=pf(q(p  ))(q{pb))(q(pe))  (c.io) 

+<5p=p[,=pe(q(p))<5p'=pc=pj:(c](p  ))(q(pa))(q(pc))  (C.ll) 

~>r^p=pa=Pd(q(p))^p'=Pb=pe{q(p  ))8p"=pc=pf{c]{p  ))  (c.i2) 

+Sp=pa=Pe(q(p))(q(Pb))  (q(pc))(q(pd))  (q(pf))  (c.i3) 

+6p=pa=Pf(q(p))(q(Pb))  (q(Pc))  ( q(Pd )>  {q(pe))  (C.14) 

+tp=pb=pd(q(p)){q(Pa))(q(Pc)){q(Pe))(q(Pf))  (C.15) 

+Sp=Pt=P/(q(p))(q(Pa))  (q(Pc))  (q(Pd))  (q(pe))  (£.16) 

+Sp=Pc=pd(q(p))  (q{pa))  {q(pb))(q(pe))  {q(pf))  (C.17) 

+Sp=pc=Pe{q(p)){q(pa))  (q(pb))(q(pd))(q(pf))  (c.is) 

+8p=pu=pb(q(p))(q(Pc)){q(Pd))(q(Pe)){q(Pf))  (C.19) 

+fip=p.=pc(q(p))(q(Pb))(q(Pd))(q(Pe))(q(Pf))  £.20) 

+tp=pi=pc(q(p))(q(Pa))(q(Pd))(q(Pe))(q(Pf))  £.21) 

+sp=Pi=pe{q(p)){q(p*)){q(Pb)){q(pc)){q(pf))  (c.22) 

+6p=Pll=Pf(q(p))(q(Pa))(q(Pb))(q(Pc))(q(Pe))  (C.23) 

+6p=pti=pf(q(p))(q(Pa))(q(Pb))(q(pc))(q(Pd))  (C.24) 

+tP=pa=pb=pc(q(p))(q(Pd))(q(pe))(q(pf))  (c.25) 

...  +  other  order-3  terms  (C.26) 

Jr5p=pa=pb=pc=pd  (q(p))(q(Pe))(q(Pf))  (C.27) 

130  ot^er  or<4er-4  terms  (C.28) 

...  +  order-6  terms  (C.29) 


The  total  number  of  partitions  (203)  in  the  summation  above  is  given  by  the  6th  Bell 
number.  Below  we  categorize  the  terms  by  their  order,  which  in  this  case  we  define  as  the 
largest  number  of  common  pixels  in  each  partition.  We  will  see  that  we  can  approximate 
the  pseudo- variance  to  reasonable  accuracy  by  including  a  particular  subset  of  the  Order-2 
terms.  Note  that  all  terms  in  the  analogous  decomposition  of  Equation  (C.3)  are  also  in 
Equation  (C.2)  and  hence  cancel  them  in  Equation  (C.l). 

C.1.2  Order-2  Terms 

These  terms  have  a  pair  of  common  pixels.  We  distinguish  between  the  following  two  cases: 

Case  1:  conjugate  pairs  These  terms  are  given  in  (C6)-(C12).  The  common  pixels  belong  to 
the  conjugate  pairs,  i.e.  p  :=  pa  =  Pd,  V  :=  Vb  =  Vo  anc^  V  :=  V?  =  Vf>  respectively,  which 
also  appear  in  the  pairwise  case.  Let  us  assume  that  the  visibilities  associated  with  a q,  <x>2, 
and  co 3  are  all  roughly  equal,  i.e.  7  =  71  =  72  =  73-  Consider  the  sum  of  the  first  of  these 
terms  (C6)  over  the  focal  plane,  i.e.: 


P=Pa=Vd,Vb,Vc,Pe,Pf 

x  (qjpe^e-^iqipfjje-^+^Pf 

=  E 

P-=Pa=Pd  Pb  Pc 

£  A p/w^  Vfe~i{-wM^f 
Pc  Pf 

(C.30) 

The  magnitude  of  this  expression  is  ( Nnpn 5 ^213)-  In  the  case  where  all  visibilities  are 
equal  to  this  expression  becomes  Nlipn5')'4.  Terms  (C7)-(C8)  follow  this  same  form. 

Now  consider  term  (C9): 
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(?(p)X?(p,))(?(Pe))«  ^ 


E 

P=Pa=Pd.P'=Pb=Pe,Pe,Pf 

x  (cj(pf))e~i(Wi+CV2)Pf 

=  E  AP  E  AfLA^'V.^/ 

P-=Pa=Pd  p'-=Pb=Pe  Pc 

x  ^Ap^-^+^Z 
Pf 

(C.31) 

The  magnitude  of  this  expression  is  N2pn4rY2,  which,  under  the  assumption  of  equal 
visibilities,  becomes  Njf)n4y2.  Terms  (ClO)-(Cll)  follow  this  same  form. 

Case  2:  mixed  non-conjugate  pairs 

We  now  consider  the  terms  in  which  a  pair  includes  two  distinct  baselines  (e.g.  (03)- 
(C18)).  The  common  pixels  in  this  case  are  such  that  one  of  the  corresponding  baselines  is 
conjugated  and  the  other  is  not.  Taking,  for  example,  the  case  (C13)  in  which  p  :=  pa  =  pe, 
we  have: 


p=pa=Pe,pb,Pc,Pd,Pf 

x  {qipcjje^+^iqjp^e-^iqipfjje-^ '+a’2^f 

P-=Pa=Pe  Pb  Pc 

x  LAPde~ia’2Pd  e~i(w\+tV2)pf  (Q32) 

Pd  Pf 

where  7a,1_a?2  is  the  visibility  of  the  fringe  of  spatial  frequency  oq  —  aq.  Note  that  the  first 
factor  in  the  Equation  above  clearly  vanishes  if  such  a  fringe  is  not  created  by  the  beam 
combiner.  Otherwise,  assuming  equal  strength  of  the  visibilities,  the  magnitude  of  this 
expression  is  given  by  (H57a»1-a;2 717273)/  or  n5ry5  in  the  case  of  equal  visibilities.  Comparing 
with  any  of  the  Case  1  terms  and  recalling  that  7  <  1,  we  see  that  the  latter  will  dominate 
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the  former  in  the  regime  considered  in  this  paper  (n  ~  1  e3,  7  <  1,  Nap  =  31). 

Case  3:  terms  in  decomposition  Equation  (C.3) 

This  Case  includes  terms  (C5),  as  well  as  (C19)-(C24).  The  common  pixels  are  such 
that  both  or  neither  of  the  corresponding  baselines  are  conjugated,  i.e.:  the  cases  pa  =  Pb, 
Pa  =  Pc,  Pb  =  Pc,  Pd  =  Pe,  Pd  =  Pf,  Pe  =  Pf-  These  terms  are  canceled  by  their  counterparts 
in  Equation  (C.2).  □ 

C.1.3  Higher  Order  Terms 

An  example  of  a  third-order  term  would  be  the  case  of  p  :=  pa  =  Pb  =  Pd,  whose  sum  is 
given  by: 


E  (?(p))^p(?(Pc))e~,'(^^)Pc(?(pe))c"i^Pe(?(p/))e"‘'(‘,,1+W2)p/  (C.33) 

P-=Pa=Pb=Pi,PcPe,Pf 

-  L  V"'E  Aprei(u,1+CV2')f’c  y,APpe  IU2Pe  LApf  e-i{cv1+w2)pf  (C.34) 

P~Pa=Pb=Pi  Pc  Pe  Pf 

The  magnitude  of  this  expression  is  given  by  h47o73-  For  equal  visibilities,  we  have  n474. 
Clearly  this  term  will  be  dominated  by  the  Case  1  terms.  Other  high-order  terms  exhibit  a 
similar  sharp  attenuation  which  allows  them  to  be  neglected  for  practical  purposes.  □ 
Given  the  relative  strengths  of  the  terms  shown  above,  we  will  retain  Case  1  terms 
(C6)-(C12),  which  after  summation  and  factoring,  yield  Equation  (4.41).  Note  that  this 
amounts  to  the  following  approximation: 


Vfizeau  (G) 


i= 1 


i= 1 


(C.35) 


Applying  analogous  analysis  to  the  Fizeau  covariance  yields  the  following  approxima¬ 
tion: 


v{Gi,Gl)fi: 


n  ( ^4)-  n  ( ^(4 ) 

keinj  keinj 


n  (zi)  n  (4)  (c.36) 

lei\J  mej\l 


Note  that  this  matches  the  corresponding  pairwise  expression  in  Equation  (4.28).  Substi- 
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tution  then  yields  the  final  Fizeau  covariance  expressions  in  Equations  (4.43)-(4.44). 

C.2  Proofs  of  Proposition  4.4.2  and  Corollary  4.4.3 

C.2.1  Proof  of  Proposition  4.4.2 

Proposition  4.4.2:  For  a  valid  RSC  array,  the  columns  Arx  and  Arv  form  a  basis  for  the 
two-dimensional  nullspace  of  Coc- 

Proof:  To  see  that  the  two  columns  mentioned  belong  to  the  nullspace  of  Coc,  note  that  each 
solution  set  to  the  noiseless  version  of  Equation  (4.49)  above  remains  valid  after  replacing 
each  Qij  with  dfj  =  9jj  —  z  ■  (r,  —  ry). 

We  then  need  to  establish  that  these  two  vectors  span  the  entire  nullspace  of  Coc  by 
showing  that  this  nullspace  is  two-dimensional.  Suppose  we  have  a  vector  w  which  is  in 
the  nullspace  Coc.  This  is  equivalent  to  the  either  of  the  following  conditions:  Mgw  =  0, 
or  Mow  E  ker(Cmc).  The  former  condition  is  not  possible  since  the  columns  spanning  the 
subspace  K  are  linearly-independent.  The  latter  condition  is  equivalent  to  the  condition  that 
M,jW  G  K  f]  L.  It  is  well-known  fact  that  for  any  two  subspaces  K  and  L,  we  have: 

dim(K  n  L)  =  dim(K)  +  dim(L )  —  dim(K  +  L )  (C.37) 

We  established  in  Section  3.5.1  that  dim(K  +  L)  =  rank( M)  =  d  +  N  —  3  for  a  valid  RSC 
system.  Also  dim(K)  =  d,  since  the  K  is  spanned  by  d  linearly-independent  columns.  Finally 
dim(L)  =  N  —  1.  Substituting  into  Equation  (C.37),  we  see  that  dim(K  n  L)  =  2.  □ 

C.2.2  Proof  of  Corollary  4.4.3 

Corollary  4.4.3:  For  a  valid  RSC  array,  the  mapping  Coc  is  injective  up  to  an  image  shift. 
Proof:  Proposition  4.4.2  showed  that  the  nullspace  is  comprised  of  linear  combinations  of 
vectors  Arr  and  Ary.  Note  that  ArY  and  Arv  are  simply  scaled  versions  of  the  x-  and  y-spatial 
frequency  vectors  in  the  array,  respectively.  Flence  adding  linear  combinations  of  these 
vectors  to  a  particular  RSC  phase  solution  merely  produces  phase  ramps  in  the  Fourier 
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domain,  which  are  equivalent  to  translations  (or  shifts)  in  the  image  domain.  In  other  words, 
the  mapping  Coc  is  invertible  up  to  an  unknown  image  shift.  □ 

C.3  Proof  of  Proposition  4.4.5 

We  begin  the  Proof  with  the  following  Lemma: 

Lemma  C.l:  The  final  2  columns  of  form  a  basis  for  the  nullspace  of  Coc. 

Proof:  This  follows  from  the  fact  Coc  is  rank-deficient  by  2,  and  standard  properties  of  the 
right  singular  vectors  comprising  V,r  in  the  SVD.  (Bretscher,  2001)  □ 

Note  that  the  error  resulting  from  application  of  the  pseudo-inverse  C+  to  the  unwrapped 
vector  of  closures  will  be  given  by: 


2nea  =  C+(2ne*h)  (C.38) 

Let  us  express  the  spatial  frequencies  measured  by  an  array  as  two-element  vectors  of 
the  form  {cox,coy).  Let  X  be  the  d  x  2  matrix  containing  these  spatial  frequencies.  Note  then 
that  the  phase-wrap  error  will  manifest  itself  merely  as  an  image  shift  if  and  only  if  this 
error  is  a  (modulo-27r)  phase  ramp,  i.e.  there  exists  a  2-element  shift  vector  z  and  an  integer 
vector  k  which  together  satisfy: 


2rcea  -  2/rXz  =  2nk  (C.39) 

Substituting  from  Equation  (C.38)  we  obtain: 

C+  (2zre)  -  2/rXz  =  2/rk  (C.40) 

Dividing  through  by  2n  we  obtain  the  equation:  C+e;*  —  Xz  =  k.  Note  that  each  element 
of  C+  can  be  expressed  as  some  rational  number  p.  Similarly  we  first  assume  X  contains 
rational  spatial  frequencies  with  greatest  common  denominator  qx.  Then  we  can  multiply 
through  by  the  least-common-multiple  (LCM)  of  the  {qt}  and  qx  to  obtain  a  system  of 
equations  whose  coefficients  are  guaranteed  to  be  integer  (i.e.,  we  have  a  linear  Diophantine 
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system).  Let  this  LCM  be  denoted  as  l.  Then  we  have,  after  rearranging  terms, 

ZXz  =  Z(C+eJ-k)  (C.41) 

We  now  wish  to  determine  conditions  under  which  there  exist  vectors  k  and  z  satisfying 
this  overdetermined  Diophantine  system.  Applying  the  Smith  Normal  Form  decomposition 
(c.f.  Theorem  2.2)  to  the  matrix  IX  this  time,  and  noting  that  rank(X)  =  2,  we  have: 


Dx  =  Ux(/X)Vx  (C.42) 

where  Ux  and  Vx  are  unimodular  matrices  of  size  d  x  d  and  2x2,  respectively,  and  Dx 
is  a  rectangular  diagonal  matrix  whose  entries  are  zero  below  row  2. 

If  we  left-multiply  Equation  (C.41)  by  Ux  on  both  sides,  we  obtain: 

ZUxXz  =  ZUx(C+e;-k)  (C.43) 

Using  Equation  (C.42)  and  the  fact  that  Vx  is  a  unimodular  (and  hence  invertible)  matrix, 
we  can  then  write: 


DXVX xz  =  l (UXC+  e*h  -  Uxk)  (C.44) 

We  are  now  in  position  to  prove  the  main  result  of  this  section,  which  is  preceded  by  the 
following  Lemma: 

Lemma  C.2:  Given  wrap-invariance,  the  (column)  vector  UxC+e;*  =  Uxetr  has  integer 
entries  below  row  2. 

Proof: 

Given  that  the  elementary  divisors  of  Coc  are  all  1,  we  know  there  exists  an  integer  vector 
ko  such  that: 


4  =  cock0 


(C.45) 


Substituting  Equations  (C.45)  and  the  pseudo-inverse  definition  in  (4.63)  into  Equation 
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(C.38),  we  obtain: 


etr  =  V^E+UjU^Vjko  (C.46) 

Noting  that  UA  is  orthogonal,  this  equation  can  be  simplified  to 

e<7  =  (Vcr  -  N)Vjk0  (C.47) 

where  N  is  a  matrix  of  the  same  size  as  \a,  which  is  zero  except  for  the  last  two  columns. 
These  last  two  columns  are  identical  to  those  of  \a,  and  hence  by  Lemma  4.1  comprise  an 
orthogonal  basis  for  the  nullspace  of  Coc.  Noting  the  orthogonality  of  Vt7,  this  can  be  further 
simplified  to: 


e„  =  k0  -  NVjko  (C.48) 

Recall  from  Proposition  4.4.2  that  the  canonical  basis  for  the  nullspace  of  Coc  is  given 
by  vectors  Arv  and  Ary,  which  denote  the  vectors  containing  the  x-  and  ^-coordinates  of 
the  baselines  in  an  array,  respectively.  Note  that  each  of  the  two  non-zero  column  vectors 
{ v/t},  k  E  1,2  in  N  can  be  expressed  (via  a  simple  change  of  basis)  as  a  linear  combinations 
of  the  elements  of  the  canonical  basis  as. 


Vfc  =  fljWi  +  fl2W2  (C.49) 

Hence  we  have: 

Uxe^  =  Uxk0  -  UxNVjko  (C.50) 

Since  Ux  is  an  integer  matrix,  the  first  term  is  clearly  integral.  Let  us  then  examine  the 
second  term,  and  in  particular,  the  product  UxN.  By  substitution  from  Equation  (C.49),  we 
see  that  for  any  of  the  two  non-zero  columns  vy  of  N,  we  have: 


UxVjt  =  fllUxWj  +  «2UxW2 


(C.51) 
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Now  note  from  Theorem  3.5.2  that  Ux  is  a  matrix  which  annihilates  all  the  spatial 
frequencies  in  the  matrix  X  below  row  2.  But  these  spatial  frequencies  are  in  fact  the 
contents  of  the  two  columns  wi  and  W2  (up  to  a  uniform  scaling  factor).  Therefore  the 
column  vector  in  Equation  (C.51)  is  zero  below  row  2.  This  means  in  turn  that  the  second 
term  in  Equation  (C.50)  is  zero  below  row  2,  and  hence  that  U x e,r  is  integral  below  row  2 
(since  the  first  term  is  integral).  □ 

To  utilize  Lemma  C.2,  we  first  re-arrange  the  Equation  (C.44)  above  so  that  it  reads: 

yDxV^z  -  UxC+ej;  =  -uxk  (C.52) 

Let  v  =  jDxV^z  -  UxC+e*.  Note  that  since  Dx  is  zero  below  row  2,  the  entries  of  v 
below  row  2  will  be  equal  to  those  of  (— UxC+e;*),  which  are  integers  by  Lemma  C.2.  Now 
consider  the  first  and  second  entries  of  v.  Let  f  be  the  vector  containing  the  fractional  parts 
of  the  first  two  elements  of  vector  UXC+  e;*,  and  let  A  be  the  invertible  matrix  consisting 
of  the  first  two  rows  of  }DxVx  1 .  Without  loss  of  generality,  choose  z*  =  A  1  f  so  that  the 
fractional  part  f  is  annihilated,  leaving  only  integer  elements  in  the  first  two  entries  of  v. 
Hence  we  now  have: 


v  =  -Uxk  (C.53) 

with  v  ensured  to  contain  only  integer  elements.  Since  Ux  is  unimodular,  the  vector 
k  =  — Ux  1  v  will  be  integral.  We  have  thus  found  a  pair  (z*,k*)  with  integer  k  which 
satisfies  the  Equation  (C.44).  Since  Equation  (C.44)  is  related  to  Equation  (C.41)  via  a 
unimodular  (and  hence  invertible)  mapping  Ux,  invariance  is  hence  proven.  □ 

C.4  Minimum  Cycle  Basis  Proofs 

Lemma  C.4.1.  (Liebchen  and  Rizzi,  2005):  Let  Li  be  the  Horton  family  of  a  graph  in  which  there  is 
a  unique  minimum  path  between  each  pair  of  nodes.  Let  C  be  a  cycle  of  G  which  is  not  in  Li.  Then 
there  exists  a  minimum  path  PU/V  between  nodes  in  u  and  v  of  C  which  is  internally  disjoint  from  C. 
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(a) 


(b) 


Figure  C.l:  Illustrating  Lemma  C.4.1  (a)  and  Lemma  C.4.2  (b) 

Proof:  Choose  an  edge  ab  in  C  connecting  nodes  a  and  b  (c.f.  Figure  C.la),  and  another 
node  u  also  in  C.  Since  C  is  not  a  member  of  'H,  either  the  path  Pua  or  Pui,  in  C  must  be 
sub-optimal  (as  otherwise  C\ab  would  be  a  shortest-path  tree  rooted  at  vertex  u).  □ 

Lemma  C.4.2.  Liebchen  and  Rizzi  (2005):  All  oriented  cycles  in  a  minimum  cycle  basis  of  a  directed 
graph  D  are  in  the  Horton  family. 

Proof:  Suppose  we  have  found  such  a  minimum  cycle  basis  B  =  C\, ...,  Cf/  Cft  and  it  was 
obtained  by  applying  the  greedy  algorithm  to  the  complete  set  of  cycles  of  D.  Furthermore 
without  loss  of  generality  assume  that  Cf  is  the  first  cycle  which  is  not  a  member  of  the 
Horton  set. 

We  know  there  exist  two  nodes  u  and  v  in  Cf  such  that  the  shortest  path  PU/V  between 
them  is  internally  disjoint  from  Cf.  Let  C\  and  Ci  denote  the  two  cycles  in  Cf  U  Pu,v 
distinct  from  Cf  and  with  opposite  orientations  on  PU/V,  as  shown  in  Figure  C.lb.  Clearly 
w(C\)  <  w(Q )  and  zv(C2)  <  io(Ct),  and  therefore  both  Q  and  C2  can  be  expressed  as  linear 
combinations  of  the  cycles  in  {Ci, ...,  Cf_i}.  Note  that  Cf  =  Ci  +  C2,  which  implies  then 
that  Cf  can  be  expressed  as  linear  combinations  of  the  cycles  in  {Ci, ...,  Cf_i}.  But  this  is 
impossible  since  the  greedy  algorithm  chose  Cf  as  linearly-independent  of  these  cycles,  and 
hence  we  have  the  necessary  contradiction.  □ 
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